texts: Can a Mosaic Dataset have rasters with different vertical units?

texts:
I am creating a Mosaic Dataset (ArcMap 10.3) that contains LiDAR DEMs from a number of different US states.

The data is in different coordinate systems (UTM, State Plane meters, State plane feet), different cell sizes (1m, 1.5m, 1.8m, 3.5 ft), and different image formats (Arc GRID, FGDB, IMAGINE, TIFF). Additionally where there isn't any LiDAR, I want to use the NED 1/3 arc second data instead (http://viewer.nationalmap.gov/basic/). The NED data is stored in Geographic coordinate system and the horizontal units are in decimal degrees. I have set the coordinate system of the Mosaic dataset to a UTM zone in the middle of my study area.

Most of the DEM data is in meters, but some of the DEMs have vertical units in feet. I want the Mosaic dataset to be consistently in meters (horizontal and vertical). Should I convert the vertical units to meters before I add the data to the mosaic dataset? Or is there something I can do within the mosaic dataset like a function that can work on a subset of the rasters.


texts: Creating shapefile with quantity of points using ArcGIS Desktop

texts:
I need a shapefile able to show the amount of houses in every 100m2 of a region. I don't want any density calculation, I need just the amount of point in each cell in a new shapefile.

How can I do it?


texts: Select by attributes appears to not work correctly

texts:
This table:

ca     ca1
12333  123330

doing this query in Select_analysis:
in the where parameter.

ca = 'substring(ca1,1,5)'

gives empty feature layer

and this:

 ca <> 'substring(ca1,1,5)'

gives a layer identical to the original:

 ca     ca1
12333  123330

Shouldn't it be the opposite since in the first case the where is True and has the 5 first numbers in ca1 similar to ca?


texts: Configuring custom pop-up window in ArcGIS Online using field hyperlinks?

texts:
I'm working on configuring a custom pop-up window in ArcGIS Online and have been using the "Create Link" feature to add hyperlinks to the pop-up window. I've been adding hyperlinks based on the name of the field in the dataset that contains a URL (e.g. {URL}).

In instances where there is no hyperlink associated with a given field (i.e. the field is empty), I've been finding that the pop-up window sometimes displays this small blue line that if clicked, will direct to a blank page. In other instances where there is no hyperlink associated with a given field, the blue line doesn't appear. In an ideal world, I wouldn't want the link to appear at all so that it is evident that there is no data to be displayed.

Does anyone know how to remove the line?


texts:
I have this code that is mostly from an ArcGIS JS API sample:

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">

    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title>Create Map and add a dynamic layer</title>
    <link rel="stylesheet" href="http://js.arcgis.com/3.13/esri/css/esri.css"/>

It works properly.
Now I copied the content of http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Specialty/ESRI_StateCityHighway_USA/MapServer in an html file and saved it in C:\inetpub\wwwroot\myweb\usa.html I started my iis service and tried to use this url http://localhost/myweb/usa.html instead of above mentioned url for ArcGISDynamicMapServiceLayer. It doen't work and it throws this error: Resource interpreted as Script but transferred with MIME type text/html: "http://localhost/myweb/test2.html? Am I missing something? This may be basic question but the answer is important for me.


texts:
This is the error I get:
GDAL command output:
ERROR 4: `file:///C:/Users/User/OneDrive/Raredirt/Home%203%20Data/home%203%202011%20jademungplanting.csv?type=csv&maxFields=10000&detectTypes=yes&xField=Longitude&yField=Latitude&crs=EPSG:4326&spatialIndex=no&subsetIndex=no&watchFile=no' does not exist in the file system, and is not recognized as a supported dataset name.
Process returned error code 1

I'm trying to create a raster of the elevation from this .csv data.

and these are my input parameters.


texts: Using MapFish to Print vector with filter based style?

texts: Is it possible to define an OpenLayers WFS style, which is based on filters, to MapFishPrint request?


texts: Determine which points fall within a 3D polygon and export this subset to a flat file

texts: I am working in ArcMap 10 and will have a 3D representation of a zone of transmission.  I would like to determine which GPS points fall within this polygon and then export this subset out to a flat file.  I have visualized info using 3D analyst but have not queried 3 D as before in the 2D world.  Thank you.


texts: Spawning a process from ArcGIS ToolValidator

texts:
I attempted to spawn a separate process from the ArcGIS ToolValidator script.  Unfortunately, it launches another instance of ArcCatalog.

What is going on behind the scenes that causes ArcCatalog to be launched, when I'm only trying to execute a second python script in a separate process and get the results?
How can I launch a separate process, wait for it to end, and get the results (standard out)?

Don't worry, I'm not trying to do heavy lifting from inside the tool dialog.  I was just trying to bypass the background processing window when calling GetCount_management.  That problem was solved, but I'm still curious why ArcCatalog gets launched.


texts: QGIS Cloud: WFS error Download of capabilities failed: Protocol "" is unknown

texts:
I published a map with several vector layers from QGIS 3.4 to QGIS Cloud (with the plugin). I activated WFS in the project properties. However, when I try to connect to WFS in QGIS desktop, the connection is established and I see the layers in the browser panel. But when I try to add them to the layers panel, I get an error message: Download of capabilities failed: Protocol "" is unknown".
How could I solve this problem?


texts: Check if function-argument is valid vector-geometry in Python GeoPandas

texts:
I wanto to build a function (let's call it clip) that accepts a vector-file (can be a shapefile, a geopackage....) and a raster. How can I make sure the argument that the users gives to the function is a type that for example gepandas can read? I though about something like (completey pseudocode):

def clip(shape, raster):
    if driver(shape) not in geopandas.drivers:
        raise TypeError
    else:
        clip(raster, shape)

I know this is kind of random pseudocode, but maybe some understand what I mean and has an idea.

When I try:

import geopandas

shape = "no_shape.sph"

try:
    shape = geopandas.read_file(shape)
except CPLE_OpenFailedError:
    print("that's not a the path to a vector geometry")
except DriverError:
    print("that's not a valid vector geometry")

I get:

CPLE_OpenFailedError                      Traceback (most recent call last)

fiona/_shim.pyx in fiona._shim.gdal_open_vector()

fiona/_err.pyx in fiona._err.exc_wrap_pointer()

CPLE_OpenFailedError: no_shape.sph: No such file or directory

During handling of the above exception, another exception occurred:

DriverError                               Traceback (most recent call last)

<ipython-input-17-7ed5ed6090fb> in <module>
      5 try:
----> 6     shape = geopandas.read_file(shape)
      7 except CPLE_OpenFailedError:

And I have no idea how to catch that exception and print a message that this is no valid shape (or any other vector supported by geopandas).


texts: Points on Map overlapping/covering city and town names

texts:
I have a map with many points, these points map incidences of landmine accidents in Colombia.  There are thousands of points and cover up parts of city and town names, in some cases the entire lables are covered.  How can I get the city and town names to be on top of the data layer, i'd rather the point be overlapped by the city name rather than the other way around.

Note, this problem is only when zoomed out, once the user zooms in on specific areas the points are small enough to not be a problem.

I am mapping on CartoDB, I imported my basemap from Mapbox Studio.


texts: Finding metadata about symbology in *.lyr file, data-source missing?

texts:
I have this lyr file which have a data source I don't know where is. When I set a new data source, it's the lyr file with all the symbology that gets replaced with the default values from the data source, no the other way around.

Is there any way to get the coloring codes (which has to point to values in a table from the data source) in the lyr file?


texts: Since upgrading to ArcGIS 10, my application doesn't rescale the map when the window is resized

texts:
I am upgrading a .NET application from ArcGIS 9.3 to ArcGIS 10.3. I've noticed that the map control behaves differently: the map no longer automatically rescales itself when the window is resized. To put it another way, when you resize the map the scale stays the same even though the extent changes. I would rather that the map rescaled when the window is resized, as it did under ArcGIS 9.3.

The commentary below from ESRI explains that this change is by design.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/What_s_new_for_map_display_and_navigation/00qp0000001w000000/

  Less redrawing in data view

  In previous versions of ArcMap, if you changed the size of the ArcMap
  display while working in data view, either by resizing the ArcMap
  window or by docking/undocking/resizing a dockable window, by default
  your map was completely redrawn to fit inside the available display
  area. So the scale changed and the extent stayed the same (although
  you may have seen some extra geographic coverage based on how well the
  extent fit inside the new shape of the display area). At version 10,
  the default drawing behavior in data view has been changed so that
  when the display size is changed, your map is no longer completely
  redrawn to fit inside the display. Instead, the scale stays the same
  and the extent will change. If you make the display bigger, you'll see
  a larger geographic extent, and vice versa.

  This has the performance advantage that the portion of the display
  unaffected by resizing doesn't need to be redrawn at all. For example,
  if you close a docked window, only the portion of the display that was
  obscured by the window needs to be redrawn. It is also easier to work
  with the display because geographic features on your map don't move
  around as you dock and undock windows. Features on your map remain in
  the same location in the display until you manually pan or zoom the
  map.

In summary: I want the ArcGIS 10 map in my .NET application to rescale when the window is resized, just as it did under ArcGIS 9.3.


texts: OpenLayers flip coordinate order for BBOX WMS request

texts:
I have a WMS layer which have a different axis-orientation to the one that OpenLayers sends out in the BBOX-request.

Openlayers uses the BBOX (minx,miny,maxx,maxy): BBOX=1438239.124213878%2C7895639.273745567%2C1448023.0638343806%2C7905423.213366069

But to get my map-tiles I have to swap the axis-order to this (miny,minx,maxy,maxy):

BBOX=7895639.273745567%2C1438239.124213878%2C7905423.213366069%2C1448023.0638343806

The x and y should change order in the request. Is there any simple solution to this?

I´ve tried to read other questions, but most requests is to change the coordinate system, which I´m not interested in.

Also, the same question I have is posted here: How change axis/coordinates/xy order for a particular layer in a Map?

But that is referring to OpenLayers 2.


texts: How can I change the size of the infowindow in CartoDB?

texts: I'm working with info windows, but I'd like them to be a bit bigger than the defaults in CartoDB. To display pictures and additional information. How can I make it so they are bigger? Is there any way in the GUI to do this? Is there a more advanced way that has some template I could easily apply?


texts: Printing subset of ArcMap layers in black and white?

texts:
I am working in ArcMap

Is it possible to print a select few layers in black and white and have one or two layers print in color?


texts:
Is it possible to get the exact Python code corresponding to a specific run of a tool in ArcToolbox?
The reason I ask is that I want to replicate the function in a Tool in Python, and adapt it to my needs.
I have tried some md tool for ncdf, which works in the GUI but not in python if I write the Python code by hand.


texts: How is it possible to publish data with ArcGIS Server without creating ArcSDE service or using ArcSDE?

texts:
I am still relatively new to the ESRI server techologies. I have server1 which runs ArcGIS server 10.0 and when I look at which processes running on this server I see ArcSOM.exe, ArcSOCMON.exe and many ArcSOC.exe (probably around 50 instances)

I have another server (server2) which runs 6 instances of giomgr.exe 11 instances of gsrvr.exe and 2 Imgrd.exe instances.

There are also 2 DBMS Servers which hold sde-geodatabase data.

My question is how have I been able to publish data from a file geodatabase stored on a partition on server2 although there is no SOC, SOM, DBMS or AGS running from it?


texts:
When I export to PDF or SVG, the layout changes and doesn't look the same as in QGIS. How do I get it right when exporting?

How it looks in qgis and how it should be:

When exporting to PDF:

And SVG:


texts: Moving an un-docked ArcToolbox from secondary screen to the main monitor

texts:
For ArcGIS 10.3.1 on MS Windows 7 (moving from multiple screens to a single screen) sometimes the ArcToolbox dialogue box remains on the non existent secondary screen. How do I move the un-docked ArcToolbox dialogue to the main monitor, without deleting the normal templates.

I've tried: Alt-Tab. Shift Space M. etc. But nothing seems to cater for ArcGIS's special little dialogue boxes. Any ideas?


texts:
I want to make an interactive web mapping application, which will grid and contour user data. My plan is to implement R scripts with Shiny, but the only R package I can find to do this is gstat. And with gstat, the only option for gridding that I can see is Kriging.

Is Minimum Curvature available in some CRAN package, or is my plan completely unrealistic and I should use something like GRASS?

Thanks for your help


texts: Exporting Data Using Attribute Table Widget In ArcGIS Web AppBuilder?

texts:
I'm using web app builder for ArcGis and in the configuration of Attribute table widget I did enable exporting data as CSV file but when I try to export it I get something like this in new page in the browser:

What I want is to download the file when the user hit export as CSV.

Is there any better way to export data from Attribute table ?


texts: Batch script converting and resampling orthophotos

texts:
I would like to convert more than 300 orthophotos produced at scale 1:8.000 with a pixel size of 0.8m X 0.8m from Ecw format to Tiff format and after resampling pixel values.

My aim is to batch convert the input format (Ecw) to output format (Tiff) and after to resample the pixel size (values) 10mX10m.

How could I right this script using gdal to batch process this raster files.

I'm using FWTools on Windows XP

Thanks


texts: Layer misalignment with UK Ordnance Survey data and OpenLayers

texts:
I am in the UK and my QGIS projects are based on UK Ordnance Survey mapping data and are set up in EPSG:27700. Whenever I have tried dropping an OpenLayers layer (Google Satellite/Bing Maps Aerial) in via the OpenLayers plugin (v0.90) it is always offset by 10 metres or so from my mapping data (rendering it useless).

I have tried turning on 'on the fly' reprojection, reprojecting in Google Mercator (EPSG:900913), starting a new QGIS project in Google Mercator projection, etc etc but the result is always the same.

I don't suppose anyone can describe a workflow that would result in a QGIS project containing an OpenLayers layer correctly aligned with (UK Ordnance Survey) mapping data?

(Apologies if this is the wrong place to post this type of query but I am inclined to think it is a case of me being stupid rather than a bug in the software.)

----Edit----

Update... The problem appears to be with QGIS, as an export of a correctly georeferenced vector map to KML as suggested by nhopton below also results in a misalignment with the satellite imagery.


texts:
I need to make a date format in QGIS 3 but, in the program I need to use,
QgsEditorWidgetSetup('Date',{'yyyymmdd'}), I want to get this:

But I don't know how to use QgsEditorWidgetSetup to make the date looks like above image. Can you help me, please?


texts: How do I generate code from ESRI GPServer wsdl?

texts:
I'm trying to use svcutil to generate code from this GPServer wsdl.

So I run this command:

svcutil GPServer.wsdl

Which generates ESRI_Currents_World_GPServer.cs.

The code has comments saying the code generator "requires additional schema information", is there some other file I can pass to svcutil that has this additional info?

[System.CodeDom.Compiler.GeneratedCodeAttribute("System.ServiceModel", "4.0.0.0")]
[System.ServiceModel.ServiceContractAttribute(Namespace="http://www.esri.com/schemas/ArcGIS/9.3", ConfigurationName="GPServerPort")]
public interface GPServerPort
{

    // CODEGEN: Parameter 'Result' requires additional schema information that cannot be captured using the parameter mode. The specific attribute is 'System.Xml.Serialization.XmlArrayAttribute'.
    [System.ServiceModel.OperationContractAttribute(Action="", ReplyAction="*")]
    [System.ServiceModel.XmlSerializerFormatAttribute()]
    [System.ServiceModel.ServiceKnownTypeAttribute(typeof(Patch))]
    [System.ServiceModel.ServiceKnownTypeAttribute(typeof(Field[]))]
    [return: System.ServiceModel.MessageParameterAttribute(Name="Result")]
    GetJobInputValuesResponse GetJobInputValues(GetJobInputValuesRequest request);

I get the same results with

svcutil GPServer.wsdl /ser:XmlSerializer

I'd like to implement the interface in a web service and have it appear to clients (like arcmap) as if it were arcgis server.  Has anyone done this?

Update:
This article addresses the issue, but the solution is to change your WSDL so that it's compatible with svcutil.  I can't do that, can I?

Update 2
It appears what I'm trying to do is referred to as "WSDL first design".  Here's a good discussion.  This description could apply to ArcGIS desktop as well as Microsoft office:

  And this it may be counter intuitive
  for some people. The client
  application, in this case, Microsoft
  Office, specifies the on-the-wire
  contract, the WSDL. Lots of people
  have a server-centric design
  perspective, and assume that the
  server defines the contract. That
  often makes sense, but not in this
  case. Because there are so many
  deployments of Microsoft Office out
  there, it makes sense for the client
  to define the contract.

Update 3 (Simplified summary)
Esri's REST API is open, yet incomplete.  Some features supported by SOAP are not supported by REST, like cancelling a GP job.  Since the wsdl is public, it should be possible to write a mapservice that implements both the SOAP and the REST API.  Does anyone know of a tool that helps me implement a mapservice that adheres to Esri's SOAP contract?


texts:
I am trying to create a hiking map. I use print composer to print a section of the map and let print composer create a UTM grid with 1000 map units distance. after i print to pdf and then to paper (A4), I measure using a ruler that the the distance between 5 grid lines is not 200mm as expected but 201mm. the features on the map is also stretched by 1 mm.

when i change the grid option in print composer to use mm instead, the grid is 200mm as expected so this is not a printer setting.

How can I fix this "stretching"?

Version: QGIS 2.18 windows and  2.18 ubuntugis
project CSR: WSG_84 / UTM32N
Map location: UTM zone 32V
scale in canvas and print composer: 1:25000


texts: Adding an empty PostGIS table as a layer in QGIS?

texts:
In QGIS 3.2 (and 3.0), when trying to load a PostGIS layer into QGIS using a Python function as follows:

def create_pg_layer(table,schema,geom_column='geom'):
    uri = QgsDataSourceUri()
    # set host name, port, database name, username and password
    uri.setConnection("localhost", "5432","","","")
    uri.setDataSource(schema, table, geom_column)
    vlayer = QgsVectorLayer(uri.uri(False), table, "postgres")
    return vlayer

It returns a None layer object, and fails to load the layer when the PostGIS table is empty, but succeed when the table contains a non-zero number of rows with non-empty geometry.

This is annoying because a PostGIS table can be treated as invalid (or valid) depending on its content.

Is there a way to allow an empty PostGIS table to be loaded into QGIS in above?


texts: How to interactively get the coordinates of a point in QGIS

texts:
I have a shapefile with a single layer of type point.

It is in the Google CRS (EPSG:900913).

When I use the 'Identify Feature' tool it gives me lots of information about the point clicked on, but I can't see the coordinates. Is there a way to find out the lat/lng of a point interactively?


texts: How to set image for Picture / IPicture for ArcGIS Desktop?

texts:
I am struggling with ESRI's Picture / IPicture Java documentation for ArcObjects 10.4.

I just need to grab the image file from my project's source (src) folder, and cast that as a Picture.
Here's what I've tried:

File pictureFile = new File("myImage.bmp");
Picture pictureDisp = (Picture) pictureFile;
//  FAIL: java.io.File cannot be cast to com.esri.arcgis.support.ms.stdole.Picture

or

Image image = ImageIO.read(getClass().getResource("myImage.bmp"));
Picture pictureDisp = (Picture) image;
// FAIL: java.awt.image.BufferedImage cannot be cast to com.esri.arcgis.support.ms.stdole.Picture

I've tried the above scenarios with IPicture as well.
I have no idea by what process on can set the Picture (or IPicture) in ArcObjects Java.

I can't use a file path, this project will be deployed as an add-in, and not all users will have access to the same file paths. I need the file embedded in my program. I just want my image file to be accessible via Picture or IPicture.

I've also tried creating via BmpPictureElement, with zero success.

http://help.arcgis.com/EN/sdk/10.0/Java_AO_ADF/api/arcobjects/com/esri/arcgis/support/ms/stdole/Picture.html

http://resources.esri.com/help/9.3/arcgisengine/java/api/arcobjects/com/esri/arcgis/support/ms/stdole/IPicture.html#selectPicture(int,%20int[],%20int[])

http://resources.esri.com/help/9.3/ArcGISEngine/java/api/arcobjects/com/esri/arcgis/carto/BmpPictureElement.html


texts: What are mystery green lines in Layout View produced during ArcObjects script execution?

texts:
Could any of ArcObjects guru tell me, what are these green lines on my Page Layout and what to do, or not to do, so they do not appear? See image attached.

I do used ArcObjects from Python. My script is converting data frame grids and graticules into graphics. As result I receive Group Element from which I need to extract single elements (text, lines etc.) in order to play with them. However, after this operation (which is carried out for 3 separate grids) I receive those mysterious green lines, which are not part of GraphicContainer. Actually, green lines appear only after conversion and extraction from one- graticule (WGS).

There is some guidance here about how to get rid of them, which works (cut and paste or delete and undelete all page elements). However, I would still like to avoid them in the first place. Also, after some configurations of cutting and pasting I figured that green lines do vanish when those 3 image elements (lower right corner) are deleted and restored.

Code suspected to be responsible for this:

pGCL.Reset() #reset graphics containere
pElem = pGCL.Next()
while pElem:
 pEP3 = CType(pElem, esriCarto.IElementProperties3)
 if pEP3.Name == '' and pEP3.Type == 'Group':
  pGE = CType(pElem, esriCarto.IGroupElement3)
  pGE_no = pGE.ElementCount
  for i in range(pGE_no):
   if pGE.ElementCount > 0:
    pGEe = pGE.Element[0]
    pGCL.MoveElementFromGroup(pGE, pGEe, 0)
   else: break
 pElem = pGCL.Next()

I know there is some discussion going on, whether GIS.SE is appropriate place for such technical, software specific questions. I apologize for this, however, I did not found satisfactory answer anywhere else.


texts: Merging unconnected line segments using ArcGIS Desktop?

texts:
I am looking for a way to connect line segments that are unconnected. It is a shapefile for a network of rivers and all the lines are slightly disconnected.

Is there something like aggregate polygons but for lines?

Edit
My river data is not continuous but I would like it to be so I can run simplification tools on it for a generalized map.

The lines don't always touch so when I try to generalize, it doesn't work. I want to connect the line together automatically.


texts: How to permanently sort attribute table records without new output?

texts:
How to permanently reclassify "Sort" the attribute  table records based on a field value ?
I try the geoprocessing tool Sort (Data Management), but this tool always need to write a new output feature or table with a correct values, what I need is to permanently sort the value on the same table or feature class.


texts: Is there a way to automatically displace colliding labels?

texts:
I have a PostgreSQL table, which i'm loading on QGIS Essen, that relates many rows to the same object and all this rows have distinct values in the field i'm labelling the objects with. Hence, when I turn on the labels and check the "allow colliding labels", my labels collide (of course!).

I want to be able to draw all labels and actually be able to read them so i'm thinking about automatically displacing them.

I've thought about having a field with a serial number for those rows that are related to the same object and then using it as the displacement "unit", but with a 50k> table, it seems like a tedious job.

Maybe there is an alternate solution?


texts: What's the status of Postgis 2.0 topology support in QGIS?

texts:
Having followed the Postgis developments over the past 6 months I have noticed a bunch of features which seem to be incredibly cool especially in the area of topologies and associating calculations on things like intersects, snapping,...

What I am trying to figure out is if anyone has or is working on any plugins which turn these features into native controls within qgis?


texts: Obtaining images to Google Drive from Google Earth Engine Python API

texts:
I have a dataframe that contains latitudes and longitudes of the centroids of rectangles that I'd like to obtain images of. These rectangles are created via the function OutputBBox, which takes in the latitude and longitude of these centroids and displaces them. This creates the respective bounding box/roi for which I would like to obtain the image for.
The code runs. The status shows "ready". However, the images do not export to my drive.
How can I export the Sentinel-2 image to the Google Drive that I want that correspond to each of the latitudes and longitudes present in the object?
I have enough space in my Google Drive folder.
def OutputBBox(lat,lon):
  north = 5/110.54 + lat
  south = lat - 5/110.54
  temp  = 111.320*np.cos(lat/(2*np.pi))
  east  = 5/temp+ lon
  west  = lon - 5/temp
  return (west,south,east,north)

for i in range(0,len(survey_buffer)):
  lat = survey_buffer.iloc[i,survey_buffer.columns.get_loc("latnum")]
  lon = survey_buffer.iloc[0,survey_buffer.columns.get_loc("longnum")]
  bound = OutputBBox(lat,lon)
  roi = ee.Geometry.BBox(bound[0],bound[1],bound[2],bound[3])
  collection      = ee.ImageCollection('COPERNICUS/S2') # searches all sentinel 2 imagery pixels...
  col_filter      = collection.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 20))
  col_filterdate  = col_filter.filterDate('2017-01-01' ,'2017-12-31')
  col_bound       = col_filterdate.filterBounds(roi)
  image           = col_bound.mosaic()

  task = ee.batch.Export.image.toDrive(
        image = image,
        region = roi,
        description = 'test1',
        folder = 'satellite-poverty-mapping',
        scale = 10
        )
  task.start()
  print(task.status())


texts: GIS and Geodatabase with Reminder Function

texts:
a friend of mine, a landscape gardener, asked me wether it is possible to represent his garden (mainly big trees and a few flower fields) in a GIS so he could administrate these on his computer by just selecting certain objects and then see a number of attributes of the object of his choice. I told him that that was no problem...

...the next think he asked me I found more interesting and usefull!

...he also wants to have reminders that tell him, when he has to water certain trees or flowers for instance or cut trees that stand close to pathways through his park.

this sounded very usefull to me but after 30 mins of google I haven't really found anything usefull in this matter.

Does anyone know if there is a way to maybe create Trigger functions in a postgresql database that triggers a reminder at certain dates or intervals?!


texts: Planet or regions offline maps between two zoom levels

texts:
I need offline maps to be included in the application I am developing. I decided to use marble https :// marble.kde.org/
to do this. Marble uses openstreetmap tileserver (among others).
What I need is to be able to access maps offline because my application hasn't got Internet access.
The first I did was to download tiles I wanted from marble.

The problem is that I don't want the openstreetmap map style, the only thing I need is coasts and relief, so I decided to create a tile server to supply tiles to marble.
I did that following these tutorial
https://ramneekpaulsingh.wordpress.com/2015/01/22/testing-your-tile-server-step-1/
http://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/

My tile server is working well because I can connect to the address http:// yourserveraddress/osm_tiles/0/0/0.png, and the gis database contains information after the osm data importation, because if I make the request
psql -d gis -c "select name from planet_osm_point where place='city';"

I got a list of names.
The renderd daemon is working fine to, because I can render manually:
./render_list -m default -a -z 0 -Z 10

But after reading the questions here, I realised that maybe it wasn't necessary, because the only thing I need is the generated tiles.

I have two questions:

Despite of the fact that my tile server is working right, because I can connect to the address http:// yourserveraddress/osm_tiles/0/0/0.png, and   I don't know how to ask it for requests. I have tried to change the downloadUrl line in openstreetmap.dgml marble configuration file for openstreetmap tile server, but my tile server seems not to be aware of the fact that marble is making requests. How can I test the rendering is being performed and it is right?
Do I really need to create a tile server? Can I render a specific area between two zoom levels and generate the tile png images in order to copy them to the host in which my application is going to be launched and marble will be able to access them?


texts:
I am doing the following
1- Go to “WEB” open “OpenLayers Plugin” and then open “Openstreetmap”
2- In the “OSM place search” I enter the name of the place I want to see on the map
3- After selecting the map, click on “create mask from selected object”
4- I go to “Raster” and cut through the mask layer
5- After completing the whole process when the recording part arrives, it gives an error ('url' is not recognized as an internal or external command, and also gives the following error “The following layers were not created correctly. C: / Users / ………………… ”)
How can I solve this problem? I just want to keep the part that I have selected on the map and for the little I know it is not letting me record it, I tried several things like running as an administrator etc.
By the way, is there any possibility on the map that I have selected to take the name of the streets and towns?


texts: Problem executing Mean Shift Segmentation: could not open the output layers

texts: I have a problem executing the Mean Shift Segmentation in QGIS with Orfeo Toolbox. I put all the parameter to execute and then, nearly 5s after appear the message that SEXTANTE could not open the output layers.


texts:
I am processing a LandSAT image for geological purposes. After some processing there are so many sandy dots in the image. To get clear geological formation as a polygon I want to remove small sandy pixels. After 'majority' filter, by 3x3, 5x5, 9x9 many of them still remaining.

If you have any better idea for how to remove them the sandy, please share!


texts:
Can I set up a variable such that the line parameters are based on this variable - for example, so I can set up line width, offset, label, and label rules based on a variable, defining the variable as a particular field of data (or an expression using multiple datasets)?

For my purposes I work with traffic projections - datasets I work with typically have speed/delay/volume for AM/PM in multiple projection years (2020/2025).

A particularly useful plot I use is a bandwidth plot (similar to the size graduation style), a visual representation whereby higher road volumes are shown with wider lines - however unlike the size graduation style, a road has two directions, and as such the lines are offset from the line position so as not to overlap the opposite direction. I also work over big areas so I don't want to show low values, so I filter these out. I may also need to use one colour for an increase in traffic on a road, and a different colour for a decrease. By this stage I have set 5 or 6 style elements based on my data, and if I want to look at the delays rather than volumes I have to manually adjust each of these style elements

I'd like to make it so I can set up my style generically, say set var="2020_AM_Volume" or var="2025_PM_Volume" or var="2025_AM_Volume"-"2015_AM_Volume", and refer to 'var' throughout my style setup - is this at all achievable?


texts: Spatial Join 3D point datasets: Ecological Marine Units (EMUs)

texts:
I have a global set of sampling points (N = 500+) that I would like to attribute with the Ecological Marine Unit (EMU) classifications near which they fall. The EMU dataset and my sampling dataset are both points datasets, and 3D, i.e. they have XYZ data.  How can I extract the EMU attributes to my sampling dataset?
I am working with the EMUGlobal package (downloaded from ESRI) in ArcGIS Pro, and ran a Spatial Join of EMU dataset to my sample point dataset. Because they are both point datasets, and 3D, I ran Match Option "Intersect 3D" with a search radius (I experimented with various radii). Every time I run the Join, all the joined fields are "0," even when my search radius is extremely large and would certainly encompass nearby EMU points.
I would like to know how to attribute my 3D point dataset with attributes from this other 3D point dataset.


texts: Display a georeferenced DEM surface in 3D matplotlib

texts:
I want to use a DEM file to generate a simulated terrain surface using matplotlib.  But I do not know how to georeference the raster coordinates to a given CRS.  Nor do I know how to express the georeferenced raster in a format suitable for use in a 3D matplotlib plot, for example as a numpy array.

Here is my python code so far:

import osgeo.gdal
dataset = osgeo.gdal.Open("MergedDEM")
gt = dataset.GetGeoTransform()


texts: Finding close road centerlines to given trajectory

texts:
I am working on trajectory prediction model for autonomous vehicle using previous trajectory and road information.

I have lane centerlines. I need to filter them. (Removing unnecessary lanes). I filter most of the unnecessary lanes by making Region of Interest(ROI).

But I need help with filtering lanes in opposite direction and perpendicular lanes in road intersections.

  My current method (using Python shapely)

def filter(candidate_cl, trajectory):
    # This function is called after lanes are filtered with ROI

    check = lambda traj_start, traj_end, line_start, line_end:
    distance(traj_start, line_start)<=distance(traj_end, line_start) and distance(traj_start, line_end)>=distance(traj_end, line_end)

    return [ line for line in candidate_cl if check(trajectory[0], trajectory[-1], line[0], line[-1])]

Basically checking if start of trajectory is close to start of line and end of the trajectory is close to end of end of line

Here are the examples

  Description: light green and red lines are trajectory, purple line is
  selected lane(green dot is start of lane and purple is the end) and green-yellowish rectangle is ROI

In first example filtering good enough, But in example 2, It doesn't remove some lanes from intersection

  Is there any way existing method to find for this problem or any tips to improve my current method?


texts: A python tool like GNU's diff for shapefiles?

texts:
I often use the GNU diff tool to compare two files or even two folders for differences. I was wondering if there is a Python based tool out there that does basically the same for shapefiles?

In the first place I am interested in attribute comparison but to find also the differences in geometry would be very handy.

Maybe someone has a script at hand I could adopt from?

I do not have access to ArcGIS so it needs to be a Python standalone solution.


texts: Use Boost Rtree with geographic coordinates

texts:
In order to find nearest nodes I'm using boost  rtree as follows:

typedef boost::geometry::model::point<double, 2,  boost::geometry::cs::geographic<boost::geometry::degree>> RtreePoint;
typedef boost::geometry::model::box<RtreePoint> RtreeBox;
typedef std::pair<RtreePoint, node_id_t> RtreePointValue;
typedef boost::geometry::index::rtree<RtreePointValue, boost::geometry::index::rstar<16>> Rtree;

// building tree
for(const auto& node: nodes) {
  rtree.insert(std::make_pair(RtreePoint { node.latitude, node.longitude }, node.id));
}

// then search for nearest neighbours for particular node
std::vector<RtreePointValue> closest;
const auto& query = boost::geometry::index::nearest(RtreePoint { node.latitude, node.longitude }, limit);
std::copy(rtree.qbegin(query), rtree.qend(), std::back_inserter(closest));

At first glance it seems to work but nodes are not sorted exactly by distance to the node. Can it be improved? In general can I use boost Rtree for geographic coordinates? If not, what can I use instead?

Boost version is 1.72.

Note: I eventually resolved it by converting nodes to cartesian coordinates.


texts:
I've been working on a plugin in python for qgis. I'm using Eclipse-pydev with remote-debug.
I believe I'm set up correctly; I've followedthe Cookbook instructions completely. When starting remote-debug, I get the blue banner with "Debugging connection activated", and Pydev Debug creates a branch named "unknown". The qgis API's path is also set in Eclipse.
From the debug I launch :

from qgis.utils import iface
layer = iface.activeLayer()
print "layername:", layer.name()

Doing this I get the following: AttributeError:'Nonetype' object has no attribute 'activeLayer'

Of course I've made sure to actually have a layer selected before running.
This means that the connection between qgis and eclipse isn't exactly working, but I've redone the instructions a few times now and pretty sure I've completed them word for word.

Here's a screencap of the debug branchings that I get:

And this is the correct debug branching as shown in pydev.org/

As you can see, the correct one branches the .py file
under 'unkown', whereas mine creates a new, probably disconnected branch.

How do I change that?


texts: google earth - viewing KML file in world map

texts:
I"m using google earth and a KML file that I have generated.  I would like to display all the KML points in a world map.  Google maps online does it but the print is poor.  Need a 300 dpi print which is a secondary issue not related to this.

How do i go about doing that in google earth or equivalent?


texts: How to work with a freshly imported vector?

texts:
I have a text file with points coordinates:

$ cat ControlPoints.txt
-10153967|8326376
-10621243|8176811
-11122894|8184333
-9657869|8379692
-10014696|8056719
-10014696|8056719

I import them with v.in.ascii:

> v.in.ascii input=ControlPoints.txt output=control_points
Scanning input for column types...
Number of columns: 2
Number of rows: 18
Importing points...
 100%
Building topology for vector map <control_points@PERMANENT>...
Registering primitives...

So far so good, but if I try to use this vector as target for some operation I get an error:

> v.what.rast map=control_points raster=soc_masked column=soc
ERROR: Database connection not defined for layer 1

A database connection must be set manually for this new vector, as alluded in the Wiki. However, creating a connection to the default database fails:

> v.db.connect map=control_points table=control_points -o
ERROR: Table <control_points> does not exist in database
       <$GISDBASE/$LOCATION_NAME/$MAPSET/sqlite/sqlite.db>

I am unable to use this vector. What should be the correct workflow to import a fully functional vector? Do I need to create a database table manually?


texts: How to handle simultaneous ArcGIS flex API GeometryService calls?

texts:
In the code below, when responses come back from the server for the b and c calls, I have the impression that the resposes are randomly given to b and c. The geometry service doesn't know which object is the owner of each response.

Do I need to declare a GemetryService in B and in C?

Class A{
    public static const SERVICE:GeometryService = new GeometryService(url);
}

Class B{
    function callService(){
        ...
        SERVICE.buffer(geometries, responseHandler);
    }

    function responseHandler(event:GeometryServiceEvent){
        ...
    }
}

Class C{
    function callService(){
        ...
        SERVICE.buffer(geometries, responseHandler);
    }

    function responseHandler(event:GeometryServiceEvent){
        ...
    }
}

Class D{
    function method(){
        var c:C = new C();
        c.callService();

        var b:B = new B();
        b.callService();
    }
}


texts: Reproducing FME Change Detector results in PostGIS?

texts:
We would like to use the function change detector from FME in PostGIS.

We had a try with the following SQL-Statements but the function does not work. There are no results in the end. The single querys for added/deleted/unchanged are working but the problem with the function we got are the two loaded layers.

Does anybody now if there is a possibility to create an change detector function in PostGIS?

The result should create, depending on the data, one to three new layers (added/deleted/unchanged).

drop table if exists tab_cd_2;
create table tab_cd_2(
gid serial PRIMARY KEY,
nume varchar,
nume1 varchar,
nume2 varchar,
geom geometry(polygon, 32633));

CREATE OR REPLACE FUNCTION cd_3() RETURNS trigger AS

$$ begin

for NEW.nume in select count(landkreise_projected_added.name_3)
From public.landkreise_projected_added
 Where not exists (select geom from landkreise_projected_added
 where landkreise_projected_added.geom = landkreise_projected_deleted.geom)
loop

select count(landkreise_projected_added.name_3)
into tab_cd_2.nume1
FROM public.landkreise_projected_added
WHERE NOT EXISTS (SELECT geom FROM landkreise_projected_deleted
 WHERE landkreise_projected_deleted.geom = landkreise_projected_added.geom);

begin
select count(landkreise_projected_added.name_3)
into tab_cd_2.nume2
FROM public.landkreise_projected_added
WHERE landkreise_projected_added.geom = landkreise_projected_deleted.geom;
end;

end loop;
RETURN tab_cd_2;
END;
$$ LANGUAGE 'plpgsql';

--  Trigger
DROP TRIGGER IF EXISTS fill_sp ON tab_cd_2;
CREATE TRIGGER fill_cd_3
BEFORE INSERT OR UPDATE ON tab_cd_2
FOR EACH row EXECUTE PROCEDURE cd_3();


texts: Configure Orfeo Toolbox to work with QGIS v2 Windows 7 from standalone installer

texts:
I have asked this question on the QGIS mailing list, but have not yet had a satisfactory result.

The QGIS manual at http://www.qgis.org/en/docs/user_manual/processing/3rdParty.html
is confusing about how OTB and QGIS work together and how to install them.

On the one hand, it indicates that the OTB files are installed as part of the standalone installation package. But then goes on to tell us what the file paths (in the Options and Configuration area) for the OSGEO4 installer.

My installation, on a Win 7 64 bit machine, was done via the standalone installer. I can see the various algorithms under Processing, but none of the OTB ones work. They all render an error that says that OTB is not configured correctly.

Because I haven't installed via the OSGEO4 installer, I can't use the file paths that are outlined on the above page.

All I can tell is that 'it should work' out of the box, but that's not my experience. Can anyone help?


texts:
I am using gstat 2.0-6 to simulate spatially autocorrelated environmental fields, but find that predict gives sometimes a fatal error that crashes R.
Below is a small subset of the data that illustrates the problem
library(gstat)
library(sp)

Longitude <- c(37.75550, 116.33333, 58.98333, 62.38000, 69.55000, 57.21667, 47.60000, 96.46000, 59.35000, 63.15000, 45.04000, 66.00000, 119.88000, 64.31000, 37.86667)
Latitude <- c(22.284800, -18.783333, 19.700000, 24.610000, 9.783333, 18.250000, -32.650000,-46.020000, 16.483333, 12.050000, 0.940000, 24.870000, -46.000000, 24.780000, -37.633333)
Temperature <- c(27.937857, 26.765201, 25.939665, 26.830605, 28.694280, 25.937617, 20.283506,  7.818690, 26.674647, 27.716699, 26.797488, 26.370757,  9.474865, 26.628669, 18.240512)

env <- data.frame(Longitude = Longitude, Latitude = Latitude, Temperature = Temperature)

coordinates(env) <- ~Longitude + Latitude
proj4string(env) <- CRS("+proj=longlat +datum=WGS84")

env1 <- env[1:10,]

v <- variogram(Temperature ~ 1, env)
fit.v <- fit.variogram(v, model = vgm(c('Sph', 'Gau', 'Mat')), fit.kappa = TRUE)
spMod <- gstat(formula = Temperature ~ 1, data = env, beta = mean(env$Temperature), dummy = TRUE, model = fit.v, nmax = 100)

v1 <- variogram(Temperature ~ 1, env1)
fit.v1 <- fit.variogram(v1, model = vgm(c('Sph', 'Gau', 'Mat')), fit.kappa = TRUE)
spMod1 <- gstat(formula = Temperature ~ 1, data = env1, beta = mean(env1$Temperature), dummy = TRUE, model = fit.v1, nmax = 100)

spMod

data:
var1 : formula = Temperature`~`1 ; data dim = 15 x 1 nmax = 100 beta = 23.47407
variograms:
        model        psill    range kappa
var1[1]   Nug    0.2623637     0.00   0.0
var1[2]   Mat 1629.3724465 12177.49   2.5

spMod1

data:
var1 : formula = Temperature`~`1 ; data dim = 10 x 1 nmax = 100 beta = 24.45988
variograms:
        model        psill    range kappa
var1[1]   Nug    0.5133904     0.00   0.0
var1[2]   Mat 1289.2962296 17142.98   1.5

predict works with the subset:
predict(spMod1, newdata = env1, nsim = 2)

[using unconditional Gaussian simulation]
             coordinates     sim1     sim2
1     (37.7555, 22.2848) 71.08162 65.14627
2  (116.3333, -18.78333) 33.52503 64.69961
3       (58.98333, 19.7) 63.94617 68.41731
4         (62.38, 24.61) 61.18118 66.69742
5      (69.55, 9.783333) 60.03236 66.89867
6      (57.21667, 18.25) 63.24242 67.85094
7         (47.6, -32.65) 63.68552 68.92589
8        (96.46, -46.02) 42.53681 64.05367
9      (59.35, 16.48333) 66.93134 66.77631
10        (63.15, 12.05) 64.10832 67.14007

but causes a fatal error with the complete dataset.
predict(spMod, newdata = env, nsim = 2)

What could be the cause (and solution)?


texts: libfreetype-6.dll lost when using GRASS command line

texts:
Update#1：

Following the instruction of Sorin Călinică, I have installed libfreetype-6.dll using the OSGeo4W Network Installer, and the command d.rast can launch successfully.

I am using GRASS 6.4.3 (installed with osGeo4w), and I have this question：

Whenever I want to use the d.rast command, I got the error message about system error of d.rast.exe that I lost libfreetype-6.dll：

How can I solve this problem? I am using GRASS 6.4.3 under windows 8.1 64-bits, please give me any suggestion!


texts: Using Landsat 5 TM 8-Day NDVI Composite?

texts:
I'm working in a project where I need the  NDVI for some areas from Uruguay for the period 1988-1993, and I found this product "Landsat 5 TM 8-Day NDVI Composite" but I'm not sure of how it was calculated because the landsat frequency is of 16 days and the NDVI value is given for every 8 days. I noticed that several Landsat products have this 8-day frequency though. so maybe this is a general thing of the Landsat products.

Here's the link of the product https://code.earthengine.google.com/dataset/LANDSAT/LT5_L1T_8DAY_NDVI


texts: Consuming real-time traffic data from Google Maps?

texts:
In order to microprocess an logistic study about the best path to take, I need to consume, somehow, real-time traffic data from Google Maps.

Is there an "open source" way to do this?

The application I'm intending to develop will start an process to consume real-time traffic data from this source, georeference the data needed for that path coordinates and then calculate, among other variables, the best path to take.


texts: Adding multiple plugins to custom menu in QGIS3

texts:
I have a set of Plugins that I have to restructure according to different sub menus in a Global Menu. For the moment, I have two plugins but each time it recreates the menu for me while I integrate the second plugin directly on the Global menu as the first plugin to Create.
class ClassGeneraleMenu:
    def __init__(self, iface):
        self.iface = iface
        # Declaration de mon Menu Global en premier lieu
        self.ambition_menu = QtWidgets.QMenu('MenuGlobal')

    # Function qui me crée des sous menus dans ma MenuGlobal
    def atrqgis_add_submenu(self, submenu):
        self.ambition_menu.addMenu(submenu)

    # Function iniGuit
    def initGui(self):
        # Insertion du Menu dans la barre menu de Qgis
        self.iface.mainWindow().menuBar().insertMenu(self.iface.firstRightStandardMenu().menuAction(),self.ambition_menu)

        # Sous Menu Par projet
        # Myprojet Sous Menu
        self.adn_menu = QtWidgets.QMenu(QCoreApplication.translate("atrqgis", "&Myprojet"))
        self.atrqgis_add_submenu(self.adn_menu)

        # CONTROLE DESIGN Submenu
        self.controle_design = self.adn_menu.addMenu('Contrôle')

        # Ajout des actions dans dans la fenetre execution
        icon = QIcon(os.path.dirname(__file__) + "/icons/import.png")
        self.animate_lines_action = QtWidgets.QAction(icon, "A.1.1 Import de données", self.iface.mainWindow())
        self.animate_lines_action.triggered.connect(self.baba_fonc)

        # Ajout de l'action dans la CONTROLE DESIGN Submenu
        self.controle_design.addAction(self.animate_lines_action)

This example is the same for all my Plugins, except that I want each Plugin to insert into an existing menu


texts: Oracle / SDE: ST_GEOMETRY or SDE Compressed Binary

texts: I read data from a SDE database (Oracle) directly. How do I know, which format is used to store the geometry data? Do I have to check the DataType for the geometry column or is this information available in one of the global tables?


texts:
I have two raster images, Boolean classified (1 or 2).  They overlap in given areas.

I need a raster that is only the overlapping sections.

Using Arcinfo 10.  Answer can only use raster methods, no converting to polygons, points, etc...


texts: Installing photo to kmz error: no module named exifread

texts: I am trying to install plugins through the Plugin Manager v2.10 and am getting the error:


texts: Too many connections on QGIS using PostGIS

texts: I'm trying to execute some queries on QGIS connected to a PostGIS. However, I'm receiving the error 'too many connections', on my company each user can only use one connection. I've searched on manual, but I did not find how to configure the maximum of connection on QGIS.


texts: How to create spatially filtered versions in ArcSDE 10.1?

texts:
I am using ArcSDE 10.1 and Oracle 11g,

I have a SDE user (parent), and two more users (ex: user1,user2) in my oracle database with required SDE permissions.
In my SDE db, i have a polygon FeatureClass (Ex: world) with a field called region (attribute values are North and South) and it is registered as versioned under sde.default.

And now i would like create, two child versions for sde.default (Ex: version1,version2) for user1 and user2 respectively and i would like to allow user1(version1) to view/edit the data related to a particular region(Ex:North), similarly for the user2(version2) to view/edit data related to the region south.

Is it possible to apply the above filter while creating versions in ArcSDE database?


texts: Using gdalUtils to reproject hdf file in R

texts:
I am now working with GEFED dataset, which is a global wild file dataset.You can download one from FTP: fuoco.geog.umd.edu; Login name: fire; Password: burnt; The file is stored into gfed4/monthly/

My task is to translate the hdf file into tif, which can be done by gdalUtils package. I think I could firstly translate it into tif file by using gdal_translate and then reproject it to wgs64 by using gdalwarp. The first step is good but I have problem in the second step. I have no idea the original srs, so gdalwarp failed.

Based on description of GEFED (can be found in gfed4 folder in that FTP), "Each product file contains seven data layers, each stored as a separate HDF4 Scientific Data Set (SDS). Units apply to each layer after multiplying by the specified scale factor. Each data layer has 720 rows and 1440 columns which correspond to the global 0.25◦ GFED grid. The center of the upper left grid cell is located at longitude 179.875◦W, 89.875◦N.". I think the original src can be written from this statement. But I have no idea how to do that and I wondering anyone can help me on this.

My code is:

library(gdalUtils)
a.path<-"G:\\CAAS_Fire\\1_GEFED\\monthly\\GFED4.0_MQ_201612_BA.hdf"
sds <- get_subdatasets(a.path)
name <- sds[1]
filename <- 'G:\\CAAS_Fire\\1_GEFED\\output\\file.tif'
gdal_translate(sds[1], dst_dataset = filename)
fileproj<-'G:\\CAAS_Fire\\1_GEFED\\output\\fileP.tif'
gdalwarp(filename,fileproj,s_srs=?????,t_srs=to.srs, overwrite = T) # I need to know s_srs.


texts: Problems with the scalebar not CRS related in QGIS

texts:
When I add one layer to the map, open it in the print composer and add a scale bar the scale bar is fine.  I can add layers (with the same CRS) and everything is fine.
Once there are two layers in my map and I open a new print composer and add a scale bar it comes out wrong (just a 0 with a line under it).  This is still the case when using the same two layers that I had in the original print composer.
Am I doing something wrong or is this a bug?
Rachel


texts: How to extract and export "raised field" LIDAR objects?

texts:
I'm working on Windows 10. In input, I have a .LAS "GROUND_only" LIDAR point cloud (4 millions points). We're interested in geomorphological microstructures called raised field (ground elevation due to cultivated areas). See picture below:

We have generated a DEM raster file (planimetric resolution 20 cms from this cloud point). Then, we have generated contour lines (10 cms equidistance). I paste below a screenshot representing these contour lines. The difference between the top of the raised fields and the bottom is often 20/40 cms:

We would like to classify our input ground point cloud in order to separate the points constituting those "raised field" and other ground points "NON ground raised field" points

Once the classification is done, we would like to be able to export the result toward a new .las or shp (ideally with useful attribute information like class_id, z value, etc.).

We would like generate a Multipoint layer only for those "raised field objects".

All the points constituting an object "raised field" must be stored into the same row (entity).

Could you tell me how to that?


texts: GeoJSON output not showing with FORMATOPTION "storage=stream"

texts:
I am trying to get MapServer to serve WFS layers as GeoJSON, using FORMATOPTION "storage=stream", but I only get an empty response in the browser.

It works if I have FORMATOPTION "STORAGE=filesystem" instead of "stream" set in the OUTPUTFORMAT block. If I test the same string (with outputformat=geojson&) with mapserv on the command line, I get the JSON output, just not in the browser.

I am running MapServer 6.4.1, and GDAL 1.11.1.

What could I be missing?

I also have GeoJSON listed in wfs_getfeature_formatlist in the METADATA section.

The outputformat block:

OUTPUTFORMAT
        NAME "geojson"
        DRIVER "OGR/GEOJSON"
        FORMATOPTION "STORAGE=stream"
        FORMATOPTION "FORM=simple"
        FORMATOPTION "LCO:COORDINATE_PRECISION=3"
        FORMATOPTION "LCO:WRITE_BBOX=YES"
        MIMETYPE "application/json"
END


texts: Measuring distace between two points with PyQGIS 2.18

texts:
Can someone help me understand what I am doing wrong in my calculation of the distance between two points?

My code is as follows:

import qgis

crs    = qgis.core.QgsCoordinateReferenceSystem()
crs.createFromSrsId(3452) # EPSG:4326
disCur = qgis.core.QgsDistanceArea()
disCur.setSourceCrs(crs)
disCur.setEllipsoid('WGS84')
disCur.setEllipsoidalMode(True)

point1 = qgis.core.QgsPoint(35.895392, -78.758094)
point2 = qgis.core.QgsPoint(35.898154, -78.680271)
disTot = disCur.measureLine(point1, point2)
print disTot, disCur.lengthUnits()

It is my expectation to get an answer of 4.36 mi or 7.02 km (from GoogleMaps) but my code gives me a distance of 0.077871997361 and a value of 3 for the lengthUnits.

What should I change in my code to the length in either meters, km or miles?
What are the units of the distance as given my code above?
Why do I get an integer value of 3 for the lengthUnits()?


texts: How to install GDAL 2.x on CentOS 7 without building from source?

texts:
I'm trying to install GDAL 2.X on CentOS7 without building from source because I'll need configure it in another servers and installing by source code will be a very tiring and time-consuming task.

I found some repositories like opengeo and ELGIS but they are for CentOS7.
On EPEL repo doesn't have GDAL 2.x available.

[root@localhost ~]# yum repo-pkgs epel list gdal*
(...)
Available Packages
gdal.x86_64               1.11.4-1.el7    @epel
gdal-libs.x86_64          1.11.4-1.el7    @epel
gdal-python.x86_64        1.11.4-1.el7    @epel
gdal-devel.x86_64         1.11.4-1.el7    @epel
gdal-doc.noarch           1.11.4-1.el7    @epel
gdal-java.x86_64          1.11.4-1.el7    @epel
gdal-javadoc.noarch       1.11.4-1.el7    @epel
gdal-perl.x86_64          1.11.4-1.el7    @epel

The version 1.11 which isn't quite the latest from 1.x major version.

Is it possible to install from any repository, or just from source?  If not, is it possible to build as a RPM to be easy installed across multiple servers?


texts: I am using ArcGIS 10.1 model builder on windows 8 to create create models for my processes. But when I open the model using Arc GIS 10.1 open model for editing with Arc GIS 10.1 running on Windows 7 it appears not to be editable. I can run the model but I cannot edit it. When I open the model in ArcGIS 10.1 running on windows 7 the model appears as on the screenshot attached below...!Does anyone have an idea what the problem might be Could? I am yet to check if I will have the same problem when I try and edit the model on another windows 8 platform... The model does not use additional tool from extensions like "spatial analyst and so on..." it only uses data management tools. Thank you in advance.


texts: How can I type an expression that give me specific text based on another field value? I added an image of example data.


texts: How to get a polyline from a polygon in ArcGIS Server 9.3?

texts:
I want to get the lines of a polygon and then get the length of each line, I'm looking on how to do this, Can anyone show me some code? Does anyone have any pointers on how to do this?

I'm using ArcGIS Server 9.3.1 and Visual Studio .NET.


texts: Change raster type gdal to raster type rasterio in Python

texts:
I have created raster using "rasterize" function with GDAL, Python
My problem is that the result is type gdal

<class 'osgeo.gdal.Dataset'>

and I would like it to be rasterio, e.g:

<class 'rasterio.io.DatasetReader'>

Is there any way I can do it without saving the gdal as TIFF and then open it again with rasterio?


texts: Adding attributes to empty shapefile from other shapefiles with different fields

texts:
I have an empty polygon shapefile with certain fields. What I want to achieve is to merge two other polygon shapefiles and then copying their attributes to the empty shapefile.

The challenge I guess is that the other two shapefiles have different field structure and naming from the empty shapefile, so this script or process need to change the names of the field accordingly.


texts: Losing information using Zonal Statistics?

texts: I produced a raster map representing estimates of the number inhabitants in pixels of 100x100 meter and I want to calculate the total of inhabitants per administrative units that are represented by polygons such as those shown in the picture. After converting those polygons into a raster layer (with a resolution of 100 meter). I'm using Zonal Statistic as Table in ArcGIS to calculate the number of inhabitants per administrative unit. The resultant table has fewer lines than the number of administrative units represented in the raster layer. Therefore I'm losing information for administrative units having an area below 28500 square meters or for administrative units with an irregular geometry such as those highlighted in the picture. The problem persists when I convert the administrative units into a raster layer with a resolution of 25 meter.


texts: How to extrude a polygon in a slope in CityEngine

texts:
I am going to extrude polygons in slope using CityEngine and would like to have horizontal roofs. As you know, when you extrude a polygon in slope, a sloped roof is generated. But in reality buildings in slope have horizontal roofs.

What is the best approach to do it?


texts: ArcGIS database feature class editor default version settings

texts:
I created an enterprise geodatabase with ArcCatalog 10.5 in my PostgreSQL database. And I created feature class in database. But when I add feature class in arcmap, I need to change the editor settings.

I uncheck the checkbox of "Edit a version of the database with the ability to undo and redo" every time. Can I set this only one time as default?


texts: Looping qgis:zonalstatistics using PyQGIS?

texts:
I'm looking for some Python + QGIS for the following:

import processing
for i in range (1,20):
    processing.runalg("qgis:zonalstatistics", "/folder/input/raster.tif",%i,
    "/folder/input/poly.shp","Var_%i",True, "/folder/output/zonalpoly_%i.shp"

The above BASH like example doesn't work, but if I substitute an integer in place of %i, the script runs correctly, so I know that the syntax is correct.

I'm a bit fuzzy on how to properly pass a variable into a string using Python and/or how to run a QGIS plugin from the built-in Python interpreter for QGIS.

Also, qgis:zonalstatistics calculates several summary statistics (mean, stddev, maximum, etc.) with no apparent means of specifying one or a few of these.

Is there a way of specifying just one, or a few?

It seems like a complete waste of processing time to calculate each if I'm only interested in one or a few measures.


texts: "Geotag and Import Photos" Paths needed for settings

texts:
Running linux ementary OS 0.4.1 Loki
        QGIS 2.18.16
trying to use "Geotag and Import Photos" and can not find the path and the .exe path  here is what I have done:

Download:
downloaded from here:https://sourceforge.net/projects/exiftool/files/

Install:
Unix Platforms
Download the Image-ExifTool distribution from the ExifTool home page
(The file you download should be named "Image-ExifTool-11.01.tar.gz".)
Unpack the distribution and make it your current directory by typing:

cd <your download directory>
gzip -dc Image-ExifTool-11.01.tar.gz | tar -xf -
cd Image-ExifTool-11.01

(At this point you may run exiftool by typing "./exiftool ".)
Test and install ExifTool by typing:

perl Makefile.PL
make test
sudo make install

Atfter adding the plugin "Geotag and Import Photos" I need to add some paths for the exiftool but can not figure out were these are:

I have installed exiftool at the home directory and it looks like this:

ANY help would be great


texts:
In Tilemill, I have a #layer that consists of a CSV with 3 columns: [Name], [Latitude], and [Longitude].

I have plotted points using the latitude and longitude, but I would like to have the Name of each point only to be visible when it is hovered over.

I am presuming that I will need to write some kind of HTML code in the Templates -> Teaser tab, but am unsure of how to go about this.


texts: Selecting Feature by SelectbyExpresion syntax in PYQGIS

texts:
I want to select features from a layer (called Road) with two fields ( A and B) with a certain condition. I want to get features with (A = MyA and B =MyB).
A and B are string.

This is my code. However, I don't know why it only selects the features according to the first part (A = MyA) and ignors the AND ??
Do you know how to use AND operator in SelectbyExpression??

MyA = T[1]
MyB =T[2]

Road.selectByExpression(" \"A\" = '{}' ".format(MyA) and " \"B\" = '{}' ".format(MyB), QgsVectorLayer.SetSelection)

selectedRoad= Road.selectedFeature()


texts: Cost distance analysis with Whitebox or other open source gis software

texts:
I need to perform a cost-distance analysis and cost allocation analysis for a tsunami evacuation model. As only open and free available data is used, it would be also good to use open source software only.

Since Grass GIS and SAGA tools do not seem to cover a cost allocation analysis, I tried Whitebox GAT (GIS Analysis --> Cost-distance-analyses).

Until now I had no success. I created a small (4x4 pixel) raster to test the tool but when I click "Run" (for the cost-accumulation) I get no result (Whitebox GAT then covers 90% of the CPU usage and seems to calculate (800mb RAM-usage), but after 15minutes I cancelled the process. I have 4gb of RAM and 4x2,1 GHz cpu).

Does anybody used Whitebox's tool for cost-distance analyses and had success?
Maybe I just don't get how the input data should be arranged, please see the attached picture - am I right? (Source cell -> all positive values)

Maybe someone knows an alternative to Whitebox.


texts: QGIS Atlas generation - Export atlas as images - only JPEG, want PNG

texts:
I can't find an option to choose image format while exporting Atlas, choices are Export atlas as images / as SVG / as PDF

The dialog that opens has no options and it always exports as not so good jpeg.

OS X 10.10, QGIS 2.6 & 2.8.1


texts: I was wondering if there is a way to angle Polygon labels in MapInfo like it's done in QGIS (see image link below).


texts: Using CURRENT map document in Python script tool?

texts:
I am attempting to use this code in a tool that I am working on:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
XMIN = df.extent.XMin
YMIN = df.extent.YMin
XMAX = df.extent.XMax
YMAX = df.extent.YMax

I then go on to do other calculations with the extent. The code works perfectly in the Python window, but because Toolbox tools are run in separate processes I am unable to use "CURRENT".

I think I could probably add a derived input variable and pass the coordinates through the validation script, but that doesn't work for what I am doing. I need to be able to do this later on in the script and add an additional temporary feature to the current map:

points = arcpy.CreateFeatureclass_management("in_memory", "NewFeature", "POINT")
arcpy.AddField_management(points, "NewField", "TEXT", "", "", 20)

I know all of my code works right, because I can paste the entire .py into the
Python window and then call Main() and run through it, but I need to share it with others in an easier way.

I would prefer the solution be compatible with 10.0 and 10.1 if at all possible, because some of the people I work with do not have 10.1 yet. I will accept a 10.1 only solution if I have to though.

What options do I have to do this?

I still can not get this to work, they have not gotten me Python Add-In Wizard yet so I am unable to test that approach. I have again tried the "Run in process" box checked and unchecked. Neither are working.

Checked gives me this: RuntimeError: CreateObject cannot open map document

Handled that exception: Object: CreateObject cannot open map document

From my testing looks like it is having the same error running in process or not.

I will still be attempting to do this within the Add-In Wizard, when they finally give me it.

The winning answer is I wasn't running it in the foreground.


texts: Is the Enhanced Vegetation Index normalized?

texts:
I am trying to figure out if the Enhanced Vegetation Index (2-band) is normalized or not? My first guess is that it is, but due to a multiplication factor of 2.4, I'm not sure if it is or not. The following paper describes the EVI2: https://www.sciencedirect.com/science/article/pii/S0034425708001971

The formula is as follows: 2.5((NIR - R)/(NIR + 2.4R + 1))

Is this an example of band normalization, such as how the NDVI index normalizes its data between -1 and 1?


texts: Using WMTS with Google maps Javascript API

texts: I am trying to learn and play around with the Google maps JavaScript API, KML, WMS, WMTS etc. I have a fair understanding of how to use KML files, Web Mapping Services, GeoJson files but I am struggling to use WMTS with the google maps API. I have gone through a lot of documentation on WMTS but could not find a working example. Could someone point me to one?


texts: GDAL: Batch create raster VRTs based on file name

texts:
I have thousands of rasters with the naming convention below.
386000_6356000(1).asc
386500_6356000(1).asc
385500_6356000(2).asc
387000_6357500(2).asc
386500_6357500(2).asc
385000_6357500(3).asc
385500_6357500(3).asc
365500_6350500(3).asc

I want to create multiple VRTs based on the suffix. For example, all files with (1) will go into one VRT, all files with (2) will go into another VRT, and so on.
I use the following code to create the VRTs but I need to be able to separate them somehow
gdalbuildvrt Raster.vrt C:\data\*.tif

Can this be done using GDAL without having to split the rasters into different folders and manually changing the code each time?


texts:
I'm trying to identify features from multiple layers. For example, I have many locations around a river (within a 100 m buffer) which I would like to identify, but instead of having to draw rectangles around river segment and identifying many unnecessary locations, I would want to only identify locations which I had previously selected (by polygon or freehand).

Is this possible in ArcMap 10.1?


texts:
I regularly have quite large mxd projects, with lots of feature layers referenced in Table of Contents (TOC). For simplicity I group them together. The problem is that I can't see the difference between a feature layer or a group (or a group of groups).

It is not only me using these documents, and it would be really useful to be able to colorize the name of groups and/or feature layers or have any other way of separate them from each other (other than having to press the +).

At the very least I'd like to have one color for groups and another for feature layers, but being able to choose freely would be a big improvement. Then I can have blue for water related stuff, red for administrative, green for... you get the picture.

Is there any such functionality in Arcmap? Or has anyone made a neat add-in that could do this?

Arcgis 10.1 on Windows 7


texts: Multiple display of vector layer in a QGIS plugin interface

texts:
def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    icon_path = ':/plugins/cropEt/icon.png'
    self.add_action(
        icon_path,
        text=self.tr(u'links crop et'),
        callback=self.run,
        parent=self.iface.mainWindow())

     # connect our select function to the canvasClicked signal
    result = self.clickTool.canvasClicked.connect(self.selectFeature)
    QMessageBox.information( self.iface.mainWindow(),"Info", "connect = %s"%str(result) )

def unload(self):
    """Removes the plugin menu item and icon from QGIS GUI."""
    for action in self.actions:
        self.iface.removePluginMenu(
            self.tr(u'&Testlink'),
            action)
        self.iface.removeToolBarIcon(action)
    # remove the toolbar
    del self.toolbar

def handleMouseDown(self, point, button):
    #QMessageBox.information( self.iface.mainWindow(),"Info", "X,Y = %s,%s" % (str(point.x()),str(point.y())) )
    self.dlg.clearTextBrowser()
    self.dlg.setTextBrowser( str(point.x()) + " , " +str(point.y()) )

def selectFeature(self, point, button):
    layer = self.iface.addVectorLayer("D:/python/cropped_area.shp", "sample_cropped area","ogr")
    selectList = []
    if layer:
        for feature in layer.getFeatures():
            # fetch geometry
                geom = feature.geometry()
                print "Feature ID %d: " % feature.id()
                layer.setSelectedFeatures(selectList)
           # show some information about the feature
            #if geom.type() == QGis.Point:
                #x = geom.asPoint()
                #print "Point: " + str(x)
            #elif geom.type() == QGis.Line:
             #   x = geom.asPolyline()
              #  print "Line: %d points" % len(x)
                #if geom.type() == QGis.Polygon:
                    #x = geom.asPolygon()
                    #numPts = 0
                #for ring in x:
                    #numPts += len(ring)
                    #print "Polygon: %d rings with %d points" % (len(x), numPts)
                #else:
                    #print "Unknown"

            # fetch attributes
                attrs = feature.attributes()
                #if attrs==row[0].strip():
                #selection.append(feature.id())
   # attrs is a list. It contains all the attribute values of this feature
                print attrs
                #layer.setSelectedFeatures(selection)
    if not layer:
        print "Layer failed to load!"

    QMessageBox.information( self.iface.mainWindow(),"Info", "Click on the feature to select" )

    #else:
           # QMessageBox.information( self.iface.mainWindow(),"Info", "No layer currently selected in TOC" )

I want to load a shapefile and want to display the properties like feature count, attribute table list in QGIS interface plugin.  Accordingly wrote the code; a portion is attached... But, on each mouse click, a new layer is getting created and no feature property is displayed even after removing the commented statements.

Can anyone suggest the error?


texts:
I created vector grids in QGIS of size 2*2km.
Now I have to add few more grids of the same size to these existing grids. Is this possible in QGIS?


texts:
I'm using the "Points to path" processing tool.  It works fine but it is not saved when I exit QGIS.  It shows up in layers panel but not on screen on restart.  I tried making the scratch layer permanent and got this error message

Could not make temporary scratch layer permanent.
Error: Creation of
data source failed (OGR error: sqlite3_open(arrows.gpkg) failed:
unable to open database file)

I don't know enough to interpret the error or correct it.  I presume I'm doing something wrong when I initiate Point to path but cannot see any obvious options to correct.


texts: Vector reprojected to the same projection in different GIS are not aligned

texts:
I have a vector file from the land use area around a river in South America. It was mapped originally in Datum SIRGAS2000 and Projection UTM, zone 23S - EPSG:31983.

As the area of study is much beyond a single UTM zone and I need to calculate areas and distances from it, I needed to change to a projection that respected areas and distances. My colleagues and I followed a set of parameters defined for Datum SIRGAS2000 and projection Albers (it does not have an EPSG code). The proj4 code for it is:

+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I reprojected it in QGIS 3 and my colleague did the same in ArcMap 10, both on Windows, but the results were different: her map was aligned with the original UTM one, but mine was displaced to the West. We already repeated the process multiple times, but with no success. Any hints on the reasons and what version of the reprojection we should trust and use?

Fig: Green: reprojected in QGIS, Black: reprojected in ArcMap.


texts: Is it a bad idea to convert a File Geodatabase table to a dictionary?

texts:
First, what am I trying to do?

I am building a Python script that should be called from an ArcMap 10.0 toolbox. This scripts uses (among other inputs) data contained in huge gdb tables. It processes the data in many different ways (hence, there are many inevitable iterations on the table).

Why would I want to convert the data?

The reason why I am converting the tables to ditionaries is that it is extremely slow to iterate through the table using SearchCursor. I have many operations that I want to do on my data and it is much simpler and faster to use existing data structures.

How am I doing it?

The usual way to use gdb tables in Python (and the one to which I have compared my method) is something along the lines of:

# Create a search cursor
rows = arcpy.SearchCursor(table)
for row in rows:
    (...)

What I'm doing instead, is that I am converting the table to dbf using TableToTable_conversion (throwing away some useless columns/rows on the way). Then I convert this DBF table to a dictionary (I was inspired by some code written by Tyrtamos, if you google "Lecture d'un fichier dbase III" you should find it). In the end, I have a dictionary containing lists of data indexed with the column names. And I'm very happy because I can do many things with it relatively quickly.

Finally,...

This method seems to be much faster. However, I am new to that game and I'm afraid to be missing something because I don't understand everything that's going on. Is it a bad idea to do that? Are there reasons not to do it other than potential memory overflow?

Please tell me if something is unclear. Thanks!

EDIT

After seeing your answers, I realize that there was one very important information missing: I'm using Arcmap 10.0.

I have run performance tests using nmpeterson's solution, modifying a few lines to make the code compatible with ArcGIS 10.0

cursor = arcpy.SearchCursor(fc, cursor_fields)
for row in cursor:
    rowList = []
    for field in cursor_fields:
        rowList.append(row.getValue(field))
    attr_dict[rowList[0]] = dict(zip(cursor_fields, rowList))

On a small table containing ~15000 rows, I measured an average elapsed time and got

make_attribute_dict -- 10.7378981873
gdb_to_dbf_to_dict-- 2.56576526461


texts: Regex for field name at ArcGIS JS API infoTemplate creation

texts:
I'm trying to make same dynamic content builder, but I'm not too experienced. This function sets the content for an infoTemplate of a sublayer(arcgisdynamicservice). It works well, if I hardcode the attribute field name like: graphic.attributes["SDE.SDE_ESRI_EP_FEATURECLASS_XX.FIELDNAME"]
I have never used regex in this type of app. Is it working at all, or did I miss something?
  function _getPopupTextContent(graphic) {
    let regex = /SDE.SDE_ESRI_EP_.*_XX.VALUE/;
    console.log(regex.test("SDE.SDE_ESRI_EP_FEATURECLASS_XX.FIELDNAME"));
    //true      running when i'm clicking on the feature
    let popupContent = "
" +
    "
" +
    "
"
      if (graphic.attributes[regex]){
        console.log("found") //not logging
        popupContent += "

AssetKKSCode	" +
                          graphic.attributes[regex]
      }

        popupContent += "
" +
                        "
" +
                        "
" +
                        "

"
    return popupContent;
  }


texts: How to union all geometry of selection in one shape with PostgreSQL?

texts: I’m using PostgreSQL with ArcGIS functions. So I trying to union all geometry of selection in one variable with sde.st_union geometry function. Then I want to insert it to another table with scene_id and shape. I use LOOP and temp variable (st_geometry type).


texts: Writing GRASS location at predefined directory

texts:
I need to write location, mapset and information in defined directory which has write permision. How can I do this? Presently the location and mapset are written at GISDBASE directory. I would like to change some other location in a directory with write permission .
Currently I have given grassrc6 as following

GISDBASE:E:ATMC
LOCATION_NAME:gcslocation
MAPSET:PERMANENT
GRASS_GUI:wxp

char* grassdbpath=getenv("GISDBASE");

after define_module();

My code like this:

G_setenv("GISDBASE",grassdbpath);
G_setenv("LOCATION_NAME","GCSLOCATION");
G_setenv("MAPSET","PERMANENT");

G_setenv("GRASS_GUI","wxpython");
G_gisinit(argv[0]);//errorm message"Location :**** not found


texts: Select query St_asBinary(rast) returns out of memory for query result

texts:
I have very large rasters (7000x3500 cells) in postgis (2.1.3) postgres (9.3.4) that I want to pull down via psycogp2 and create numpy arrays and ultimately a geotiff.

Initially I found I had to use ST_AsBinary(rast) instead of St_AsGDALRaster(rast) in order to prevent python from running out of memory.  I'm using windows 32bit python 2.7.7.

The script I wrote does exactly this:

Query the raster as binary
Decode the binary and create numpy array
Manipulate the numpy array
Save array to geotiff.

I used it occasionally over the last two months without any issues and yesterday I started getting DatabaseError: out of memory for query result

At first I thought this was a python memory issue so I tried the query in PgAdmin III and got the same result. out of memory for query result

Here is my query:

SELECT ST_AsBinary(rast), timezone, validdatetime
from my_raster_table
where validdatetime = '2015-10-1'
and timezone = -6  -- use UTC offset
and model = 'the_only_model_we_currently_have'
and parameter = 'the_only_param_we_currently_have'

After some troubleshooting one of our developers found he was able to pull the same raster down ok.  The difference between our queries was he didn't include the model and parameter in the where clause.  They are currently not necessary so I removed them and found I'm able to execute a single select OK.

SELECT ST_AsBinary(rast), timezone, validdatetime
from my_raster_table
where validdatetime = '2015-10-1'
and timezone = -6  -- use UTC offset

HOWEVER, when I kick off my script to request multiple dates I'm only able to pull down the first couple rasters before I run into the error again.

I should note that each time I query the database I create a new connection and cursor which are closed immediately after the cur.fetchall(). The script is doing this in serial.

I do not know which memory this error is referring to.  It has that feel of a memory leak, but I have no evidence to suggest that it is one. All signs on the server suggest postgres is not out of memory.  All signs on my client machine suggest I'm not out of memory.  So where is this error occurring?

TLDR: When querying a large raster via ST_AsBinary, why does 'out of memory for query result' occur and what can I do about it?


texts:
Is there a way to recover a non functioning aprx file?

In ArcMap you could use mxd doctor, but from what I've seen ESRI don't have anything like that for ArcGIS Pro.


texts: How to apply a default datum transformation in an expression in QGIS

texts:
I'm using the expression:

x(transform($geometry, 'EPSG:32632', 'EPSG:3003'))

to add coordinates fields according to different CRS in a point layer. To obtain better accuracy i need to apply a Datum transformation (epsg:1662) in the CRS conversion. I have set the option "ask for default datum trasformation if several are avalaible" in Settings > Options that it should allow to select the datum transformation but using the expression in the field calculator this function does not seem to be available (the window does not appear).

How can I overcome this problem if possible or it would be better to create a custom CRS (how?) appplying the datum transformation to insert in the expression?


texts: Updating rows based on number and position in table using ArcPy cursor?

texts:
I'm trying to update rows using the update cursor, but based on the number of rows in a sequence of three list items.  Looking at the list below, it needs to update first four rows to 'B', then the next 9 rows to 'D' and then the next 14 rows to 'A'.

I've gotten as far as the code below, but all I'm getting is an update of the first four to 'B' and then I get 5 calculated to 'D' and then another 5 to 'A'.  I think it's going back to the start of the table somehow.

RowCount = 1
lenClusters = [4,9,14] # number of rows to update in that sequence.
rows = arcpy.UpdateCursor(fcPhotos, "", "", "CapturedStatus")

for row in rows:
    if rowCount > 0 and rowCount <= lenClusters[0]: # See if the current rowCount is from 1 to 4
        row.CapturedStatus = "B"
    elif rowCount > lenClusters[0] and rowCount <= lenClusters[1]: # See if the current rowCount is from 5 to 9
        row.CapturedStatus = "D"
    elif rowCount > lenClusters[1] and rowCount <= lenClusters[2]:  # See if the current rowCount is from 10 to 27
        row.CapturedStatus = "A"
    rows.updateRow(row)
    rowCount += 1

I also can't rely on the row ID as it might have missing IDs.


texts: Avoid duplicate labels on river sections and choose placements

texts:
I am trying to label rivers's names.
The rivers include multiple  sections. I have many repeats like "fleuve Loire ".
How to avoid labeling repetition? "Merge connected lines to avoid duplicate labels" does not seem to be sufficient.

How to place the labels at the "middle" of each river (sum of sections /2)?


texts: PostgreSQL, PostGIS and GeoServer publishing 1.7 million records

texts:
I have a polygon geodata in postgresql and I am publishing out to geoserver the 1.7 million records.  How long is too long when processing?

EDIT:  Forgot to mention that I successfully published a database spatial view subset of this same data.

Side question, What is postgresql really doing during this process???

@Devdatta Tengshe, I have a polygon dataset with about 1.7 million records and I am just trying to publish the data to geoserver.  I published a spatial view from a small table and a join to the large polygon dataset.  This layer seems fine in openlayers.  I just published another spatial view from a different small table joined to the large polygon dataset and it published but will not finish loading in openlayers.


texts:
I am trying to import a WFS layer into QGIS. I tried the normal procedure, Add WFS Layer, and it does not display any boxes. When I click on Add, it throws an error:

  Layer is not valid: The layer restrictToRequestBBOX='1'
  srsname='EPSG:4326' typename='bu:Building'
  url='http://ovc.catastro.meh.es/INSPIRE/wfsBU.aspx' version='auto'
  table="" sql= nis not valid and cannot be added to the map

I also tried the WFS 2.0 Plugin and in this case, I was able to download some vectorial data which seemed to be ok, but only if I was zooming very close to extract some features. When I zoom out to get the extent I need (or even a bit more), I get the following error:

  OWS Exception returned from the WFS: Area of extension out of limits.

However, I know there is data beyond my bounding box.

Am I doing something wrong? I also changed the EPSG several times and it did not work either.


texts: Fitting many rectangular shaped polygons into larger irregular shaped polygon in PostGIS?

texts:
I am looking to try and find the optimal fit for a number of rectangular shaped polygons into a larger irregular shaped polygon and I haven't found much information about how this could be achieved in PostGIS. I have found a solution using ArcGIS Desktop at Fitting known size polygons into irregular polygons using ArcGIS Desktop? However, I don't have a licence.

Doing some research it looks like I will need to use some kind of bin packing algorithm, however most of these algorithms deal with fitting small regular rectangles or squares into a large rectangle or square.


texts: How to cluster 2 datasets in Google Maps API V3?

texts:
I'd like to produce a map that shows thousands of Markers.
I'm using Markersclusterer to show all these markers.
Yesterday, I received another huge data set (completely different) and I'd like to see its markers with the ones that belong to the previous data set (using MarkerClusterer.
So, I can have different icons for those two data sets.
Can you help me do this?
thanks


texts: Create Buffer in meters Spatialite query

texts:
I am new to Spatialite. I have following query:

select A.*
from linka as A, pointa as B
where Contains(Buffer(B.Geometry, 100), A.Geometry)

Here,

  A.Geometry = LINESTRING(23.2161 69.95237, 23.21581 69.95279)
  B.Geometry = POINT(23.21817 69.95177)

I want to create 100 meters buffer and get to know which are the link's are contained by point geometry.

I am having coordinates of WGS 84 projection system from navigation data.

The passing 100 over here, it takes as degree and creating buffer like that, but I actually want to specify the buffer distance as meters. (I don't want to pass the distance in degree value because on each earth pole degree's buffer cover not same range area. This is the big problem of geodecy)

Is there any way to satisfy my requirement in Spatialite ?


texts: How can you overlay raster data on a stamenmap in R?

texts:
I'm trying to overlay an image output from Google Earth Engine (.tif) on a basemap obtained using get_stamenmap {ggmap}.
The object returned by get_stamenmap doesn't have a coordinate reference system (calling crs(some_stamenmap) returns 'NA'). I infer that the object returned by get_stamenmap uses Google's PseudoMercator projection, because whatever it's using plays fine with feature data output from GEE
(Aside to explain what I mean but that: if I export a featureCollection from GEE, read into R, extract the long/lat data from GEE's 'javascript dictionary' format ( GEE outputs a csv with a single .geo column containing data of the form {"type":"Point","coordinates":[107.35873495689798,15.07126214132908]} ), and then throw this data on the map using ggmap(basement)+geom_point() it displays correctly.)
The raster I am trying to overlay on a stamenmap basemap has the following properties:
class      : RasterLayer
dimensions : 1800, 4008, 7214400  (nrow, ncol, ncell)
resolution : 0.08983153, 0.08983153  (x, y)
extent     : -180.0044, 180.0403, -76.77901, 84.91774  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

I have tried reprojecting this to Google PseudoMercator (EPSG:3857) with no luck. How do I fix what appears to be a CRS problem, and get my raster data on my stamenmap, when the stamenmap basemap has no CRS?


texts: Limiting Input Parameter to only File Geodatbase in Python Toolbox?

texts:
I've been working with ESRI Python toolboxes for a while but I'm still having difficulties setting an input parameter as a "DEWorkspace" and using the parameter filter to make it so only file geodatabases are visible to the end user when entering parameters.

I'm using 10.5 and this is the ESRI documentation

    param0 = arcpy.Parameter(
    displayName="Input Workspace",
    name="in_workspace",
    datatype="DEWorkspace",
    parameterType="Required",
    direction="Input")

# Set the filter to accept only local (personal or file) geodatabases
    param0.filter.list = ["Local Database"]
    params = [param0]
    return params

The problem I'm encountering is that this filter doesn't seem to actually limit what the end user can see or enter as parameter. For example when they navigate to provide a file geodatabse the dialog allows them to see and enter .lyr files.

I'm assuming that you should be able to limit what the end user can see and enter when it comes to workspaces?


texts: How can I improve the resolution of a vector grid made in north pole stereographic projection?

texts: I've been struggling with this for a while and wondered if someone could lend a hand.  I have a map which already has sampling points from centroids in a 35x35km grid that was created in north pole stereographic.  I need the grid to represent 5x5km, essentially dividing the grid by 7; however, I can't for the life of me figure it out.  I've tried to buffer, I've tried to duplicate and move the affline (large area, so it gets distorted), etc.  Densifying the lines creates equal spacing along the edges of the polygons but I can't figure out how to get them to join in lines, because of the offset of the lines due to the projection warp.  Attached is a picture of what it looks like - any ideas?


texts:
I work for a mid-sized school district and we unfortunately have to close one of our school buildings.  That means that the neighborhoods who've been sending their children to the school that is closing will need to be redistributed among the remaining schools.

Currently, I have several kml files loaded into Google Earth, including one for each grade level of students attending the school to be closed, and one for each of the remaining schools containing a polygon representing it's boundary.  Visually, this is very informative, but it's not giving us what we really need.  Hopefully, there is another tool or feature out there that somebody can point me to.

I'm hoping to find a tool or website where I can manipulate the polygons representing the other schools' boundaries, pulling them into the neighborhoods of the closed school, and have it tell me how many placemarks or pins are then inside of each polygon.  I need to experiment with different configurations to try to evenly distribute both the overall number of students, but also the number of students at each grade level.

Does such a tool or web site exist? Unfortunately, I do not have the programming know-how or time to design something myself. Alternatively, perhaps somebody has another suggestion on how to accomplish this goal?


texts:
I put this tag in a text block of the mainPage:

- !text
    text: 'Page: ${pageNum}/${pageTot}'

But it appears only the current page number!

Someone tell me knows how to solve this problem?


texts: ArcMap shows error 000539 for NDSI calculation

texts:
I want to run Raster Calculator to calculate NDSI from atmospherically corrected Landsat images, but it shows error 000539 and some other errors.
I downloaded bulk data (Landsat collection 2 level-2) from Earth Explorer and used the NDSI formula in Raster Calculator.
The calculation works on the image downloaded from EOS Landviewer but I wanted to use an atmospherically corrected image. I have included a screenshot of the output message.


texts: Simplifying a multipolygon into one polygon respecting its outer boundaries

texts:
I have a multipolygon containing a one large polygon, itself containing a number of smaller polygons, which constitute 'holes' in the largest polygon.

My goal is to simplify the multipolygon into one polygon that closely respects the existing perimeter but removes all the smaller polygons (holes)

I have been using a combination of ST_ExteriorRing and ST_Dump to get rid of the points and lines, however I can't make this work when there are multiple polygons.

The closest working solution I've found is to use ST_ConcaveHull(geometry, 0.99), however this does not follow the perimeter of the shape closely enough (it smoothes over the recesses).

I've pasted the text representation of the multipolygon below
https://gist.github.com/glennpjones/51acef849825470386f24a6c3259295d

I'm using PostGIS 2.5.

How would you approach this?


texts: Plotting new offset points in ArcGIS for Desktop?

texts:
I have a task where I have plotted positions along an alignment (see red dots in images) and I need to plot two offset points (from the chainage locations) at varying distance parallel to the alignment (see image 2). I am not sure how to go about it.

I have thought about:

working out the direction of the line and from that,
then working out what the offsets x,y should be using trigonometry,
create new points for the offset using the results of the trigonometry.

I was also thinking of:

creating a 1m buffer around the points,
using the buffer to clip the line into segments, then
copy the line parallel on the left and right side using the specified distances
then use feature vertices to points tool with mid as the point type to get the offset positions.

If there were only a few points, I would use the second method above, but unfortunately, I have hundreds.

I am familiar with Arcgis 10.1 and have a sound understanding of Arcpy and python.


texts: I would like to make a map of sinkholes. I have a  DTM ARC/INFO ASCII GRID of Asiago Plateau and and I'm using grass gis. Does anyone know what is the correct form in grass gis?


texts: Getting user's QGIS version using PyQGIS?

texts: I am trying to write a Qgis plugin, which should work on Qgis2 and Qgis3. How can I find out which Qgis version the user is using, so the right imports are made?