texts: ArcMap doesn't recognize projected coordinate system in geotiff

texts:
I am working on a NASA grant to produce gridded satellite passive microwave
data products.  Our products are formatted in netCDF-CF and we are producing
geotiffs of data layers using gdal_translate.  The data are projected to one
of three EPSG ProjectedCRS codes (6931, 6932, 6933).
The gdal developers have been very responsive with changes to make sure they
handle our data correctly.

I have been testing the resulting geotiffs with ArcMap 10.4.  I expect that
I will be able to open legitimate geotiffs and reproject them using ArcMap.
However, while the current geotiffs that I am producing seem to have
metadata that looks entirely accurate and correct, ArcMap is throwing the
following error when I try to open them.

  The following data sources you added are missing spatial reference
  information.  This data can be drawn in ArcMap but cannot be projected.

Sure enough, the data image is rendered, but the the frame says Current
Coordinate System is "No coordinate system".

The gdalinfo output for this file looks like I would expect:

$ gdalinfo NSIDC-0630-EASE2_N25km-F13_SSMI-1997061-19H-E-GRD-CSU-v1.0.nc.TB.tif
Driver: GTiff/GeoTIFF
Files: NSIDC-0630-EASE2_N25km-F13_SSMI-1997061-19H-E-GRD-CSU-v1.0.nc.TB.tif
Size is 720, 720
Coordinate System is:
PROJCS["WGS 84 / NSIDC EASE-Grid 2.0 North",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",90],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","6931"]]
Origin = (-9000000.000000000000000,9000000.000000000000000)
Pixel Size = (25000.000000000000000,-25000.000000000000000)

I would be obliged if someone can tell me why ArcMap is unable to read our data.

I can say that I worked with EPSG in 2013 to define these 3 ProjectedCRS entries, so I know that they are relatively recent additions to the EPSG data base.  And I wonder if perhaps ArcMap is using EPSG and if it needs to be told to use the most recent EPSG release?  But I'm only guessing on this, I don't have direct evidence or enough knowledge about how ArcMap is actually working.


texts: Embed a (saved) Leaflet map in an Ioslides presentation

texts:
I'm building some fairly complex maps with the leaflet package and I'd like to embed them in an rmarkdown document (specifically an ioslides_presentation .Rmd file).

Right now the code for the maps is in an .R script and it is subject to change as I tinker around with map aesthetics. I know I can cut and paste each map's code into the .Rmd presentation doc, but I want to avoid have to manually update the .Rmd file every time a change is made. I figure that if I can save the leaflet output and embed it in the presentation that will be much smoother.

I know that htmlwidgets::saveWidget allows me to save leaflet maps. But how can I embed these outputs in the .Rmd presentation?


texts: Calculate The Minimum Distance from Polygon Centroid to Edge

texts: I am looking to calculate the minimum distance from the Center of a polygon (A crescent shaped polygon would have a center inside the polygon) to its edge. I have parcel data and I am trying to calculate in my model I created how tall I can build a cell tower on that parcel. The maximum height that I can build the tower is equal to the minimum distance from the center of the parcel to its edge. Is their a simple way to calculate this? Once I have the minimum distance calculated I can complete my query for my model.


texts:
How can I design a point such as in this image below ?
How can I compose the color in RGB or hex or so on to get this point?


texts: Converting a multiband raster to a single band raster wipes off its values

texts:
I have a multi-band raster and I am using

gdal_translate -b 1 -scale 30 30 landscape.lcp band_1.tif

to get a single band in tif format. Now if I check values of band 1 in landscape.lcp, those lie between 0 -270. But if I check values in band1.tif, the values lie between -32768- 32767. Any pointer how to keep the original values in the band_1.tif file.


texts: How to access a secured ArcGIS REST API using Python and extract the data in the form of Shapefiles or JSON

texts:
I have a secured ArcGIS REST API which requires username and password to access it. I am able to access a public service and extract the data in the form of shapefile. Below is my code,

import urllib2
import arcpy
import os

myurl="http://services.thelist.tas.gov.au/arcgis/rest/services/Public"
myservice="/TopographyAndRelief/MapServer"
myparams="/8/query?where=objectid%3E50+and+objectid%3C2001&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&f=pjson"
myRequest = myUrl + myService + myParams
myRequest = myurl + myservice + myparams
response = urllib2.urlopen(myRequest)
myJSON = response.read()
foo = open("C:\\Users\\venkatac\\Documents\\jsonOutput1.json", "wb")
foo.write(myJSON)
foo.close()
arcpy.JSONToFeatures_conversion("C:\\Users\\venkatac\\Documents\\jsonOutput1.json","C:\\Users\\venkatac\\Documents\\finalShapfile1.shp")

But I need to access a secured service using the credentials and extract the data.

Can anyone give me suggestions or any idea?


texts: QGIS (maptool) Plugin Crashes on some computers

texts:
I developed a custom Plugin which inherits from QGIS MapTool.The plugin has user interfaces and can add features to the map based on a GUI. The Plugin works Fine on some computer(when they installed QGIS with OSGeo4W). But on standalone version does not work and crashes and gives the following error.

The plugin loads without any error even show the GUI but when I click on the map to place feature it shows crashed.

The problem is I can not put the codes here as they are about 9000 line codes.
Can anyone explain what is  QObject::Connect and QHashData::allocateNode

As no one understood the problem I decided to add more information:
I think there might be somthing wrong with the way I connect the signals:
here is my code for connecting signals:

 def __form_action_handler(self):
    """This Will handle the buttons action/ and field change  """

    self.main_dlg.OperationalMode_btn.clicked.connect(self._operation_mode)
    self.main_dlg.UGStructureType_cbo.currentIndexChanged.connect(self._form_info_change_handler)
    self.main_dlg.AddToList_tbtn.clicked.connect(lambda:self.add_to_structure_to_add_list())
    #self.feature_list_dlg.remove_tbtn.clicked.connect(lambda:self.remove_from_to_be_added_list_table())
    self.feature_list_dlg.remove_tbtn.clicked.connect(self.remove_from_to_be_added_list_table)
    self.main_dlg.Hierarchy_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','Hierarchy_cbo', "hierarchy"))
    self.main_dlg.Size_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','Size_cbo', "size"))
    self.main_dlg.Material_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','Material_cbo', "material"))
    self.main_dlg.ObjectStatus_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','ObjectStatus_cbo', "object_status"))
    self.main_dlg.PlantOwner_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','PlantOwner_cbo', "owner"))
    self.main_dlg.JointUse_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','JointUse_cbo', "joint_use"))
    self.main_dlg.SpaceOccupancy_cbo.currentIndexChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','SpaceOccupancy_cbo', "space_occupancy"))

    self.main_dlg.NodeHousingID_txt.textChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','NodeHousingID_txt', "id"))
    self.main_dlg.AdditionalID_txt.textChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','AdditionalID_txt', "additional_id"))
    self.main_dlg.Lat_txt.textChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','Lat_txt', "Lat"))
    self.main_dlg.Long_txt.textChanged.connect(lambda:self.__single_field_vlaue_change_handler('self.main_dlg','Long_txt', "Long"))


texts: How to get an ocean shapefile from a very detailed land area shapefile?

texts:
I am looking for a possibility to "invert" a shapefile. What I want to do is to convert a shapefile with the continents and islands into a shapefile that has water polygons. The existing ocean shapefile from Natural Earth is too inaccurate.

As an active OSM contributor, I thought using the existing shapefiles from OSM would be fine for me. These are much more accurate but very buggy from a GIS perspective. However, mapserver can deal with it so I thought these should be fine. These shapefiles, some of them you can find here
http://metro.teczno.com/#coastline
are very large.

I have created a shapefile with one polygon with world bounds (-180 to 180, -90 to 90). Then I have tried to clip out the areas from my world bounds shapefile where land areas are in the other shapefile. I tried it with QGIS, vector tools -> difference function.
After 12 hours with 100% CPU usage, I stopped it. It always displayed a progress of 100%. RAM was not an issue.

I am looking for another possibility to clip out the land areas from a "world polygon", preferably Linux command line.


texts: Getting a new column with distance to the nearest feature in R

texts:
Having two sf data frames (a is points and b is polygons),
I can get columns from the nearest feature by for instance:

c <- st_join(a, b, st_nearest_feature),

But how is the best way to get an additional column with the distance as well?


texts: Create an ArcCatalog OLE DB Connection using Python

texts:
Is it possible to create an OLE DB connection (*.odc) to a database in ArcGIS 10 , but using Python instead of the ArcCatalog GUI wizard?  I don't need to necessarily create a file on disk, just make the connection to the database.  I am trying to automate a procedure that many users may use, and do not want to rely on every user remembering to create the OLE DB connection manually through ArcCatalog.

The Python code provided by @blord-castillo in a similar question shows How to create an SDE Connection using Python, so I am looking for something similar that creates an OLE DB connection.  I cannot seem to find any reference to any Python code for this task.

Eventually, the python code sample using Make Query Table is basically what I want to do, but it assumes you already have your *.odc connection built.

Anyone have suggestions or pointers for places to look for sample code?  I am working against a SQL Server database, in case that matters.


texts:
The situation is:
There are many photos (between 900 to 3000), geotagged and I can see them on MapInfo. The photos are of certain objects and there may be 2 to 8 photos for each object. These objects are also visible in MapInfo as a separate layer. Something I do with them is I go through each one to find the best photo and delete the rest so only one photo per object exists.

The way each photo is currently being assigned to an object is by physically moving the icon over the object. The point in doing that is to ensure that the assigned photo is closest in distance to the assigned object. An algorithm later updates the tables of each layer based on a shortest distance method.

There is no field in the photo layer to say which object it belongs to, and there is no real obvious geometrical pattern between the photos and objects. The only way to tell which object the photo belongs to is by visually inspecting the photo and then doing some kind of data entry for that photo.

Here is a sample of what the map looks like.

The problem is, sometimes I make an error (perhaps by rushing or losing concentration) by thinking that two photos of the same object but from different angles are photos of two separate objects. So the consequence is that there are >=2 photos of the same object when there is only meant to be <=1 photo per object.

The desired outcome is that I can identify photos that belong to a particular object (by somehow grouping them) in a quick manner for the purpose of running a check to ensure only one photo per object exists after I have gone through and found the best photo icons for each object icon in MapInfo.

In other words, I want to group the photos some how that will make this statement true about a particular group of photos: 'this group of photos belong to one object'. So, when I've finished finding the best photo per object, I can run some kind of script to check that <2 photos exist in the table of the same object.

Manually inputting the data into a table in Excel would look something like the image below. This would work for what I want the data for but it is going to be too slow to be practical.

The question is:
How do I select multiple photos in MapInfo, which will then make them belong together and so I can later check to make sure only 1 photo per group exists after I have gone through and selected the best photos?

Please feel free to provide a comment if you have another idea about achieving the desired outcome.


texts:
We are using CartoDB as our production mapping solution.

Just wondering does it offers support for Retina screens and if so how hard can be to toggle on?


texts:
How do you add a wmts layer from the USGS basemaps found here as a mapbox-gl-js layer?:

https://viewer.nationalmap.gov/services/


texts: Does QGIS have anything like Geodatabase Versioning?

texts: I am just wondering if QGIS has anything like Esri's geodatabase versioning?


texts: How can i add custom menu in QGIS Using python

texts:
I am tried to add a custom main menu in QGIS Desktop Application. Here I am already added custom sub menu in QGIS but I want to add a custom main menu in QGIS using Python.
See the blow image.


texts: Shapely (geos) crashes during unary_union

texts:
When attempting to union multiple polygons with Shapely 1.5.6 in Python 2.7.9, using GEOS 3.4.2, I get the following error in the underlying C++ code:

Assertion failed: (!"should never be reached"), function itemsTree, file AbstractSTRtree.cpp, line 371

The polygons are in a Python list and then I call shapely.ops.unary_union(polygon_list).  The polygon list has 1998 items, all of which are proper polygons.  I have tried running the unary_union operation in a separate python session with a couple manually created polygons and it works.

Could the number of polygons be a problem, or the import order in the file (I have seen a mention of shapely needing to be imported before Basemap, but I am not using that package)?


texts: Getting a rotated shape when I change projection using GeoPandas

texts:
I have a GeoPandas dataframe that I can properly plot
f,ax=plt.subplots()
myGdf.boundary.plot(ax=ax)

If I change projection I get, a rotated geometry of 90 degree.
myGdf1 = myGdf.to_crs({'init': 'epsg:3174'})
f,ax=plt.subplots()
myGdf1.boundary.plot(ax=ax)

Am I doing something wrong?


texts: Is it possible to display a gray color GeoTIFF in ERDAS as it simply occurs in QGIS?

texts:
I have several GeoTIFFs (gray color) to be imported in ERDAS, however after googleing I could not find a way to display that color gradient as it occurs in QGIS.
Below the GeoTIFF displayed in QGIS:

And below the imported GeoTIFF as IMG and displayed as "Pseudo color" in ERDAS:

If that can be displayed correctly in ERDAS, then how?


texts: Will there be a Python API for ArcGIS Server?

texts: It seems like an obvious next step given the big push towards Python in ArcGIS Desktop at 10.0.  Are there any limitations of Python compared to javascript, flex, etc that would prevent ESRI from developing a Python API?


texts: Clipping vector polygons by previous vector polygons in same layer in QGIS

texts:
I have extremely large datasets as vector polygons in single layers.  The polygons are in an order.  I want the previous ordered polygons to clip the later polygons.  Therefore, once an area has been used it can't get used again.

Example:

polygon 1 would never be clipped
polygon 2 would be clipped by polygon 1
polygon 3 would be clipped by polygon 1 & 2
polygon 4 would be clipped by polygon 1 - 3
...
polygon 5000 would be clipped by polygon 1 - 4999 if no unused space
is left for polygon 5000 or any other polygon it gets eliminated.


texts: In ArcGIS Desktop 10 I have an arc that has violated the topology rule that the road feature should not extent outside of my boundary.  I tried to use the clip tool in editor but cannot complete the clip. I should be able to click on the boundary, then use the trim tool to clip off the extended line.  I selected the boundary, its lit up and shows in the TOC as selected, but nothing happens with the clip tool and I can't do anything with the offending road feature that is over the boundary line.


texts:
From an SRTM DEM I have calculated the slope of the terrain with values that range from 0 to 77 (degrees). I have also calculated the curvature of the terrain and the statistical characteristics of the outcome are as shown below:

What does "e" mean? Is the result of curvature reasonable?


texts: Exporting selected features to shapefile PyQGIS

texts:
I am trying to develop a script to export selected features on a layer to enable some automation.

So far I have used this formula to import a layer and run a query. This ends up with features being selected on the layer.

See result:

All good so far.

I now try to use the below script, found in a previous similar question:
QGIS Export Shapefile using PyQgis

It worked for the users in the post, however this returns the following error for me:

I am quite new to Python so am not quite sure what's wrong. It is saying name 'i' is not defined, but I'm not sure what 'i' is trying to refer to?

I'm using QGIS 2.18.20.


texts: Working with both ArcGIS and Oracle Spatial, a self-intersecting polygon renders a geometry invalid.  Why are self-intersecting polygons a problem?


texts: how to map/drape polygon layers over raster in QGIS?

texts:
I want to "map" or "drape" a 2D vector (polygon) over a raster. Obtaining a 3D vector with z-values according to the raster it's draped over.

Most of my internet findings point towards v.drape.
However, this tool seems to only work for polylines, not polygons. Furthermore I am looking for functionality similar to that of ArcGIS 3D-Analyst ("interpolate shape"), where points are inserted in every new raster cell with an according height.
Is there a similar tool or combination of tools to achieve this?

Unfortunately, converting the polygon to polylines, draping it and converting it back to polygon does not seem to work either, as the Z-values get lost.

QGIS 2.18.16


texts: OSGeo4W error - The procedure entry point ... NCSEcw.dll

texts:
There are a few posts with a similar error but none with this particular dll. When I try to run a gdal executable from the command line or try to open the OSGeo4W shell I get the following error...

In the terminal I get the error...

ERROR 1: Can't load requested DLL: C:\OSGeo4W\bin\gdalplugins\gdal_ECW_JP2ECW.dll
127: The specified procedure could not be found.

There are no duplicate dlls in the OSGeo4W folder. I have also deleted the OSGeo4W folder and reinstalled. Still get the same error. It also appears that Pathfinder Office has the same dll, but it's in a different directory. The system worked great until a couple of weeks ago with both Pathfinder and OSGeo4W installed. I'm not sure what has changed since.


texts:
I'm using GeoExt.Popup to create a popup when a user click on an WFS (point) on the map.
When the user click on a point on the map I want to send the GID (primary key from DB) to a servlet or JSP to send back HTML with all information about that point (using iframe inside the GeoExt Popup).

I don't know how to get the GID value for the point that user clicked.
Any one knows how to get an attribute from the map?

I'm using this code:

    if (popup) {
    popup.destroy();
    popup = new GeoExt.Popup({
    title: "Info:",
    width: 600,
    height: 350,
    margin: 1,
    autoScroll: true,
    maximizable: true,
    collapsible: true,
    map: map,
    location: map.getLonLatFromPixel(e.xy),
    //html: e.text
    html: "<iframe src='response.jsp?codVertices=1' width='100%' height='100' frameborder='0'></iframe>"});
    popup.show();

In the "HTML" option I want to get the GID and send it to the servlet or JSP.
It's working fine when I manually insert the GID in the html code...

Thank you very much.
Pedro


texts: Advantages of using the export acknowledgement message?

texts:
Could anyone tell me what is the advantages of using the “export acknowledgement message” in replica?

I couldn’t figure out what is the real meaning of using this tool. What is the effect of it if we don’t use it after replica in disconnected environment?


texts:
I will simplify my question to be better understood. Let's suppose my study area is consisted by settlements A and B. They both have the same population (1000 people each). Each settlement was converted to raster and then to point. Let's assume that each settlement has 100 points (each point corresponds to 10 people therefore a population field was added and filled up). The maximum radius of the walking distance is 3000 meters. Location 2 is accessible by many points from A and many from B (e.g. 100), location 1 is accessible by many from A and few from B (e.g. 60) and location 3 is accessible by none from A and many from B (e.g. 40). I want to show a raster output (visually similar to a density map) where each cell is assigned a value that corresponds to the number of points that can access it in a radius of 3000 meters.

In reality the shape of the settlements is complex and the population is different. The question is "accessibility" and I initially thought I wanted a "point density" output but I realized it doesn't show what I want. Are you aware of any tool that would help in this regard?


texts: Performing Batch Geoprocessing in ArcGIS Pro?

texts:
I have been looking around in order to add multiple files instead of 1 at a time. I currently have about 55 files that I need to keep as separate layers.

ArcGIS 10.x had the option where you can right click a toolbox and choose Batch, but that's not the case (as far as I know) in ArcGIS pro.

Is there a way for me to use the tool for multiple files in ArcGIS pro, or it's only done through python or ModelBuilder?


texts: How to vary the infoWindow size for different feature layers

texts:
I have multiple feature layers on the map, and I wish to display a different sized infoWindow for each layer.

As far as I can tell, the map has one infoWindow, which can be resized using:

map.infoWindow.resize(width,height)

Each feature layer then provides the infoTemplate which contains the Title and Content, but no sizing information (right?).

How can I show a large infoWindow when clicking on FeatureLayer1, and a small infoWindow when clicking on FeatureLayer2?


texts: Mapping ratio function over ImageCollection in Google Earth Engine

texts:
I am trying to map a function over an ImageCollection to receive the ratio in a polygon, how it compares to a certain point (pt) for each image. So all pixels within the polygon divided by the value at point pt.
I keep on receiving the error message "Error in map, Image.constant: Invalid Image.constant type." and cannot find a way to solve it.
I tried different codes, e.g. this one:
var maskL8 = function(image) {
  var qa = image.select('BQA');

  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image.updateMask(mask);
}

var start = '2020-07-01';
var finish = '2020-07-30';
var pt = ee.Geometry.Point([-49.31582,69.56833]);

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate (start,finish)
.filterBounds(pt)
.filter("WRS_ROW <= 122")
.map(maskL8);

var myB8 = l8.select("B8");

//receiving the values at pt for each image

var getB8 = function(image) {

  // Reducing region and getting value
  var value_B8 = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), pt)
    .get('B8');

  return value_B8;
};

var count = myB8.size();

var B8_list = l8.toList(count).map(getB8);

print("B8 list", B8_list);

//Ratio function

function  calculateRatio (image) {
  var ratio1= image.clip(polygon)
  .divide(ee.Number(image.reduceRegion(ee.Reducer.first(), pt))
    );
    return ratio1;
}

var RatioALL = myB8.map(calculateRatio);
print (RatioALL);

So I would like to receive the ratio for each image so I can afterwards get the standard deviation and the mean for each clipped image as output. Maybe it would also make more sense to first do the function and clip afterwards but both showed the error for me.
I included the code for getting the values at pt for each image in case that helps for the function.


texts: Exporting contour lines to x-y-z raster in ArcGIS Desktop?

texts:
I am using ArcGIS 10.5 and have a shapefile containing contour lines. I would like to get the height values at equidistant positions. In other words, I would like to interpolate and sample the contour lines. The result should be a table containing height values for the range of x-y values of the shapefile, depending on the sampling distance.

The goal is to export the data as a simple CSV file to use it in programs such as MATLAB or Python.

Unfortunately, I am stuck finding the correct tool or process in ArcGIS for this.

In the ArcToolbox I found Conversion Tools -> To Raster, but this seems to raster the contour lines themself.

Also, I found the help describing how to export a complete XY coordinate values to a text file from ArcMap, but here the XY coordinates are not equidistant and I am unable to export to TXT as the menu only allows exporting to a file with an extension of .dbf.


texts: Attribute errors with selectbylocation algorithm

texts:
Could anybody tell me what is wrong with execution of

processing.run('qgis:selectbylocation', layer2, layer, u'intersect', 0)

layer2 and layer is QgsVectorLayer object.

This is error message I got


texts: Script for batch-clipping or masking rasters in R

texts:
I need a script that can do following batch process for all files in a folder:

a) Clip or mask each raster layer (LAI product from MODIS MCD15A3) with a vector layer (vegetation types obtained from the MODIS land Use product)

b) Calculate the average for the pixels in the clipped raster

c) Divide pixel value by 10 (the MCD15A3 is LAI x 10)

d) Create timeseries with table values (date, LAI)

e) Export to CSV

f) Delete the intermediary files

I am aware that this is multiple questions in one, but i have browsed through the Q&A's here and many other places and I have been unable to find anything that will solve the task.


texts: How to transform RGBI (4 band) to RGB (3 band)?

texts:
I got 4 band RGBIs and I need to transform them to 3 band in order to be able to use them in Agisoft Metashape.
I didn't find a way to delete the 4th band (NIR), so I created 3 new images (each containing 1 band) by using the raster calculator. Afterwards I merged them. As result I got a image with 3 bands, but the first band always gets "gray" in brackets. And when I export them the outcome is a white picture.


texts: Concatenating Fields where Values Match in other Field using ArcGIS Field Calculator?

texts:
Is there a way in ArcGIS Field Calculator to concatenate fields from a column if values match in another column?

I have attached an image for reference.


texts:
When I try to open my workspace file, renamed 'FILE' here, it gives this error:

(FILE.wor:201) Found [)] while searching for [A string constant].
(FILE.wor:204) set is an invalid indentifier.
(FILE.wor:204) Unable to create map layer.
(FILE.wor:204) Unable to create map group layer.

I have tried opening the file in notepad to see if I could find a 201 or 204 but there are none, so it must be embedded somewhere in the file.
Who can tell me how to fix this?

Thanks!
Bertus


texts: Adding PNG/SVG symbols on the map on GeoServer

texts:
I am trying to add all symbols (cafe, hospital, school, and etc...) on my map on GeoServer.

For the test I added this code (to see the cafe symbol):

name: Test
title: Test
feature-styles:
- name: name
  rules:
  - filter: ${type = 'cafe'}
  symbolizers:
  - point:
      symbols:
      - external:
          url: file:/Downloads/symbols/Cafe.p.16.png
          format: image/png

But I got empty screen, even there weren't buildings, lacks... Only white screen.

Then, I tried:

name: Test
title: Test
feature-styles:
- name: name
  rules:
  - filter: ${type = 'cafe'}
  symbolizers:
  - point:
      symbols:
      - external:
          url: file:/Downloads/geoserver-2.11.1/resources/120px-Cafe-16.svg
          format: image/svg

But result is the same.

Polygon, line, road styles are SLD files, but point is YSLD.
If it's important, the order of styles is:

Polygon
Roads
Line
Point


texts: Geographically Weighted Regression with overlapping polygons

texts:
Is it statistically correct (or possible at all) to calculate a spatial regression with (partly) overlapping spatial areas?

I have market areas as the spatial unit with different sociodemographic, economic and other (firm-endogenous) variables and I´m examining their influence on the economic success (demand). To examine necessary regression assumptions a firstly need to calculate a simple OLS-regression to check for spatial autocorrelation (Moran I. test), homoscedasticity, normal distribution of standardized error terms etc. …Only from the visual perspective, when I look at the spatial distribution of my regression variable values I suppose a strong spatial autocorrelation so I probably need to implement spatial lag variables in my model. I´m using R for my calculations. To calculate the spatial weighting matrix of the neighbours I need to read the SHP-files to T. But how the spatial weighting matrix can be calculated when the areas are overlapping?

Aggregating my data to a higher level to avoid the overlapping effect is in my case not possible (too few cases for the regression analysis).
I´ve only found some papers dealing with overlapping observations in time models:

Improved Inference and Estimation in Regression
With Overlapping Observations

The Overlapping Data Problem

but I haven´t found anything dealing with spatial overlapping.

Do you know some methods, useful transformations etc. to solve this problem with spatial overlapping or is it possible at all?


texts: Creating map or converting to GeoTIFF, Bufr format?

texts:
I'm trying to work with spatial data saved in WMO BUFR format (.buf extension). I need to cut data to my extent, create map and calculate spatial statistics.

I'm aware that two tools are able to open BUFR files:

Metview (but only RedHat-based packages are available and compiling it from the source is real pain)
Grads (but to open it, one needs to prepare ctl description file for each bufr),
but none of them (?) is able to export BUFR to GeoTIFF/other GIS format.


texts: How to combine several rasters into one attribute table with no value changes in ArcGIS?

texts: I have a set of 1000 rasters that I need to copy into excel to process some data. I want to know if there is a way to combine these rasters so I have one attribute table that keeps the properties of each individual rasters. therefore when I copy the attribute table into excel I can difference the data from each raster. The rasters don't overlap with each other.


texts: How to select by attribute and loop through values to GetCount per month with ArcPy?

texts:
I would like to write a script where I select a layer by the [District] attribute and loop through the field [PickUpDate] to get the count per district per month.  My month format is m/d/yyyy, i.e. 1/2/2015.
Below I select the district in the Table View, I could easily use GetCount to get the total for my selection, but I need the count per month per District

Example:

alayer = 'PickUpColumns_12_1_2015_dups1'
hbclause = "[District] = 'HB'"
arcpy.SelectLayerByAttribute_management(alayer,"NEW_SELECTION",hbclause)


texts: I need to write a simple app for IPad that displays maps I was given. All maps are in one folder and they are shapefiles. I am able to open them in ArcGIS for desktop, however I am unable to attach them to iOS application. I've gone through the tutorial with iOS SDK, however they use .tpk files. I converted my maps into .tpk, however my iOS app crashes when I try to open that converted tpk file, unlike other tpk files provided in tutorials. Could someone give me the direction with this problem?


texts: How can I merge rasters using GDAL inside Python interpreter or script?

texts: I am in need of merging rasters within Python script or interpreter, not outside using a command line tool such as gdal_merge.py. Is there a simple way? Or is there a way of executing the command line tool from inside a python script?


texts: Pixels, distance to line (negative left, positive right) r

texts:
I am using r to do some spatial analysis.

Here is my project, I am using raster files (that contain information on land use) and a polyline (which is a border between two regions). I am analyzing land use change for each pixel according to the position of that pixel regarding the border. In a nutshell I want to do a spatial regression discontinuity analysis.

Here a visual figure:

I have been able to collect the data of each of my (17 rasters) using that code :

r1 <- raster("~path/r1.tif")
r2 <- raster("~path/r2.tif")
...
r17 <- raster("~path/r17.tif")

rasStack = stack(r1,r2,...,r17)
points <- rasterToPoints(rasStack,spatial=TRUE)

write.csv(points, file = "~path/points.csv")

This provides me with a csv file that contain 20 variables (point_id, lat, long, r1_value, r2_value, ..., r17_value).

I do not want to aggregate the pixels into bigger ones. Hence I have about 22 million pixels.

Here is my question: I am trying to figure out how to add to my csv file one variable with the distance between each pixel and the line. I want to find a way to get a negative distance if the pixel is on the left, and a positive distance if the pixel is on the right.

I found one alternative way to do it on QGIS, by (1) converting the border shapefile into a chain using chainage, (2) create a distance matrix for each pixel to the nearest chain, (3) create a one side buffer around the line border, (4) clip the pixels to this one side buffer, (5) merge clipped pixels to the previous one and set them as either being on the right or on the left of the border.

However I am trying to do it in r to avoid doing things manually, and because I am sure there must be one easier solution.
I found some functions like get.knn but did not get the expected results.

Here an image of the line and the whole raster:

Here an image of how the line intersects with the raster:


texts: Formatting numbers coming from dataset to include commas for thousands in Carto map legend?

texts:
I'm using carto.js v4 to build my map. The legend I've put together displays some of the data as numbers next to value scaled markers
As example I have:

    document.getElementById('gasDistmax').innerHTML = `

${metadata.getMax()}
(MCF)

The ${metadata.getMax()} gives me 585457. I would like to present that number as 585,457. It would be preferable to not change the actual column format in the db or add an extra column formatted as a string.

Is there a javascript way to do this, either within the carto library or one of the js.min libraries?

I've been building this from various online examples but have yet to find anything to reformat the numbers.


texts: Google maps distance based on the zoom level

texts:
I have used google maps for displaying some locations using markers, based on the user's current location. User's current location is the center point of the map. The default zoom level is 12. In default, I displayed markers within 300 miles(In my case) distance only.

In that how can I do the following,

1) If decrease the zoom level, I need to increase the distance that should display markers within the increased distance also. how much distance should I increase in miles (For e.g. Zoom level is 11) ?

2) If increase the zoom level, I need to  decrease the distance that should display markers without the increased distance. how much distance should I decrease in miles (For e.g. Zoom level is 13) ?

Please suggest me there is any built-in option to do that ...


texts: What is the best way to display a large shapefile (with polygons) as a layer on a leaflet map?

texts:
I have a leaflet map of the UK with around ~2000 points total points and a bit of metadata loading for each point with two side-panels showing the metadata. These points are drawn from several tables stored in a PostGIS database on my own VPS. Now I’d like to add layers to display some standard demographics, i.e. census, index of multiple deprivation, etc. underneath the points. The trouble is, most of the spatial demographics for the UK are associated with census output areas (OAs) and the shapefile which includes all the UK OAs (for which I'm using the CDRC's excellent aggregated shapefile) is very large, around 250MB. This is only a technical problem, as I prefer maximal granularity for these polygons so OAs are good for this. But for this to work on the web, I need to find a way of reducing the amount of data being drawn from the server.

I’ve already tried generalising this shapefile in a series of ways (using postgis, qgis, and GRASS v.generalize - which are all some variation on Ramer–Douglas–Peucker method) but even pushing this to a level of generalisation that becomes unrecognisable gets me to around a 7mb toposjon file, which is way too big to work for the web. I expect that the same user map look at this map both at a National level (zoom level 9 or so) and on a local level (zoom level 13+) so need proper fidelity for the polygons but also need to be careful not to download hundreds of MB worth of unsimplified polygons when working with max zoom.

Going back to the drawing board, from what I can tell, there are three options:

Stick with geojson, but split my shapefile into a series of smaller geo/topojson files (but will need to also have a simplified national level geojson shapefile)
Import the .shp files into my PostGIS database and then have leaflet calculate bounding box for current leaflet map instance displayed onscreen using map.getBounds(); then draw using PostGIS ST_MakeEnvelope() to filter polygons returned by the query. Possibly also run the PostGIS SQL query through ST_Simplify() if working with zoom level < 11 (or so)
Set up a tileserver to serve pre-generated vector tiles (against using bounding box) proxied from my PostGIS database using something like Mapnik,  T-Rex or Tilestache or load everything into GeoServer and serve using WMS.

I’ve scored the web and stackexchange and though most Q&A and tutorials are oriented towards working with a small geojson (usually a simplified map of the USA from Mike Bostock’s fantastic tutorials) there are a few overlapping but slightly different queries here:

Visualizing large datasets of polygons with Leaflet
Millions of polygons: displaying layers in Leaflet based on location being viewed?
and more indirectly here: What is the most common way of displaying geodata from Postgis on Leaflet?

So what’s my best (i.e. most efficient) option? Is it one of the three above, or perhaps something else? The main concern here is to be sure that this leaflet map can load efficiently (without noticeable delay) across a range of zoom levels. I'm also a tad hesitant to implement a fully-featured tiling server unless this is the only feasible option. I'm hoping someone with experience setting up a similar solution can speak to the actual optimisation trade-offs (server overheads vs. client-side bandwidth) implicated in their preferred option.


texts: Set Relative Hyperlink to open windows explorer using Actions in QGIS

texts:
I am wondering how to solve this problem I am having. I am trying to make an action in QGIS when I click on my feature, windows explorer will display the folder for that feature. I have the folder name in the attribute column, and an absolute path to the top level of the folder in the Actions.

I have looked up how to hyperlink in QGIS and have looked on the forums, but I have not been able to find a good solution to my questions. I have tried the eVis plugin and it is not a solution I am looking for. I want specifically windows explorer to open when using action on a feature, which is set using relative path.

I hope the screen capture shows up, but if not...
In Layer Properties - Feature|Actions, the Type is set to Windows, the name is Folder, the action is recorded as:
explorer C:\Users\Me\Desktop\TN\Test\Picture[%"Hyperlink"%]

How would I get the recorded action set to a relative path instead of an absolute path?


texts: Existing C++ libraries that convert ESRI shapefile to DXF/DWG

texts:
Is there any C++ library that can convert ESRI Shapefile to DXF/DWG file?

I seek a link/or sample on how this can be done with the said C++ library.

I'm aware that one can do this with ArcGIS Desktop application itself, but this is not what I want. I want to write code that can do this with my own UI.


texts: How to create a plugin GUI for QGIS 2.8.2?

texts:
I create a pluging with the pluging builder in Qgis 2.8.2. I went to the folder of my pluging : C:\OSGeo4W64\apps\qgis\python\plugins\yu in order to visualize al the files that the pluging builder created. After i launch the OSGeo4w batch file, in order to create the .ui file of my pluging by using the command:

pyuic4 -o ....

Here's the error that i have:


texts: FileGDB feature layer processing using OGR2OGR sqlite scalecoords

texts:
This is an extended question of the following link :  Multiply or Divide certain numbers on a projection file(?) of Japan

It works great with shapefiles but somehow it's not working well with FileGDB feature classes

Some of my layers are over 2GB; therefore, i had to store them as FileGDB feature class. I have installed FileGDB driver (using the latest version of GDAL) and modify the code shown in below

ogr2ogr -f "filegdb" -a_srs espg:4301 -dialect sqlite -sql "select scalecoords(geometry, 2.7777777777777777777777777777778e-7), ATTRIBUTES1, ATTRIBUTES2, from D:\NAME_VS_E.gdb\NAME_VS_E" D:\NAME_VS_E.gdb\test D:\NAME_VS_E.gdb\NAME_VS_E

and it returns an error
"Unable to open datasource'D:\NAME_VS_E.gdb\NAME_VS_E' with the following drivers.

I have installed the driver though.
It works great with

ogrinfo -so D:\NAME_VS_E.gdb

Mysterious....


texts: Linear model in R - Multiplication Expression

texts:
I have 3 numerical variables A, B and C. I am trying to create a linear model capable of predicting A. The expression that I am using is the product of B*C in order to predict A; however, when looking at the output I am not able to get my equation because I get and extra variable that I don't know what is it.

Here is my code

MyData<-read.csv("...", header = T)
head(MyData,6)
str(MyData)

#Linear Model
#Expersion A= B*C
Model1<-lm(MyData$A~MyData$B*MyData$C)
summary(Model1)

Output of str(MyData)

> str(MyData)
'data.frame':   6 obs. of  3 variables:
 $ A: num  2.5 3.4 2.7 3.6 2.5 2.1
 $ B: num  0.01 0.02 0.015 0.017 0.018 0.01
 $ C: num  0.1 0.2 0.27 0.19 0.17 0.16

Output of summary(Model1)

Call:
lm(formula = MyData$A ~ MyData$B * MyData$C)

Residuals:
       1        2        3        4        5        6
-0.03945 -0.08386 -0.13925  0.67703 -0.40055 -0.01393

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
(Intercept)          5.473      5.774   0.948    0.443
MyData$B          -222.431    454.508  -0.489    0.673
MyData$C           -26.482     36.222  -0.731    0.541
MyData$B:MyData$C 1938.961   2679.207   0.724    0.544

Residual standard error: 0.5688 on 2 degrees of freedom
Multiple R-squared:  0.6149,    Adjusted R-squared:  0.03723
F-statistic: 1.064 on 3 and 2 DF,  p-value: 0.5178

lm uses the Wilkinson-Rogers notation so "*" is an iteration, based on the output, right? is this true, how do I create my model using the product of my two variables?


texts: Esri enterprise GDB features showing as non-spatial table in QGIS

texts:
Why are Esri's Enterprise geodatabase 10.4.1 (in MS SQL Server 2012) feature classes showing as  non-spatial table/spatial type is 'no Geometry' in QGIS?

The storage properties of the feature class in ArcCatalog shows as ST-Geometry spatial type(SRID 300001).


texts: How to do Network Kernel Density Estimation Analysis in ArcGIS for Desktop?

texts:
I am looking at doing some analysis of road traffic collisions. Obviously collisions occur on stretches of roads and I was researching the best methodology to take into account the network element of it. The usual hot-spotting (Kernel Density) smooths things out too much, and I found a few articles on Network Kernel Density Estimation which is supposed to be more representative as it calculates the distance between events via the network rather than using euclidean distances. We are using ArcGIS Basic/Standard 10 at the moment. We do not have Spatial Analyst licence but we have Crime Analyst and Network Analyst.

Anybody has done this before in ArcGIS? I've googled it I can't find any ready made tools - so I was wondering if anybody had a script I could use?  Or any idea on how I could do it in a model?


texts: GeoServer HTTP header proxy authentication

texts:
I have an instance of GeoServer running on Tomcat configured so that it can be accessed through both SSL on 443 and with non-SSL on port 80. The SSL connection requires client authentication and the non-SSL allows only anonymous read access, so the user must provide a client certificate in order to gain admin access. Here is a rough diagram showing the configuration

Because of new security requirements I have to do the authentication using oracle webgate on an httpd server configured to be a reverse-proxy for tomcat. In this configuration GeoServer would be doing authentication using "HTTP Header Proxy Authentication" Although this configuration works for the admin access, I would be forced to close off the non-secure connection to the world outside of the firewall since it's easy to bypass security by adding the right property to the HTTP header (i.e. set the oam_remote_user=admin in the header and you're logged in as an administrator). Here is a diagram showing this:

The problem is that we have a user base which depends on anonymous, read-only access to consume data. I feel this should be common enough an issue that other people have dealt with this, but haven't been successful at finding any posts on this. I have dug through GeoServer documentation and tried to fumble around with the settings without luck. I would hope that there is a setting where proxy authentication is only allowed through a particular IP address or port or a range of IP/Ports, but haven't found any such setting. Is there are way to restrict it to an ip address? If not is there another way to accomplish the same thing?


texts: How to calculate with time values in QGIS?

texts:
I have two string fields in a vector layer, "From" and "To", with values in the format of 10:16:30 and 10:18:00 e.g.

I would like to calculate the difference in minutes as a decimal value, in this case 1.5, in the field calculator and write it to a real field for further calculations.

As I understood so far, I have to convert my strings to QT time type by using totime("From") and totime("To").

But now I already reached the end of my wisdom. No matter what I try, the field calculator doesn't accept any further combination of other functions (Expression is invalid).

Can anyone please tell me how to perform this calculation and receive a real field with a real minutes value, or can point me to some resources with working examples?


texts:
I have a shapefile with polygons. And I have a csv file with points.
Now I want to know for each point in which polygon it is situated. Every polygon has a number, so I want to create a file where you can find for each point the number of the polygon it is located in.


texts: Openlayers: Update scale from resolution in cached features?

texts:
I have a large number of cached vector features which I am caching for speed. However, as they need to change size depending on view resolution, I have to do various updates on them with each increment, which sort of defeats the purpose of caching.

Is there a way for cached styles to update their scale (and other options, like offSet) when resolution changes?

Here's the basic code:

const imageVectorStyle = function(image, resolution) {
  // Caching both the style and the icon, fyi.
  let style = styleCache[image.getProperties().name];
  if (style) {
    //The style exists in the cache, but since the resolution has changed
    //since caching, I have to update the scale.
    style.getImage().setScale(1 / resolution);
    return style
  }
  else  {
    const imageProperties = image.getProperties();
    let imageIcon = iconcache[imageProperties.src];
    if (!imageIcon) {
      imageIcon = new Icon({
        offset: imageProperties.offset,
        size: imageProperties.size,
        stroke: new Stroke({color: '#dee2e6', width: 1 }),
        // Here scale is set based on resolution when the style is first created.
        scale: 1 / resolution,
        crossOrigin: 'anonymous',
        src: imageProperties.src
      });
    }
  }
}


texts:
I need to color each polygon based on its area and length. I can loop through all polygons, find out their area .geometry().area() and if it is larger than 100 meters squared, than I want it to be green, otherwise red.
Is it possible? I found solution for one attribute. (Setting color of feature depending on attributes with PyQGIS?)


texts: Obtaining least number of grids (polygons) covering a separate polygon while the grids are overlapping

texts:
I have grids as polygons where the grids are overlapping. I want list of grids covering a separate polygon. The image below demonstrate the problem. The red lines are the edges of the grids and the black polygon is polygon of our interest.

Since the grids are overlapping simple intersect will generate a list where some grids add to spatial redundancy. I want to optimize size this list by removing the grids which only contribute to the spatial redundancy. For example, in the following image the skyblue polygon will intersect with four surrounding grids, but the expected output is a single grid (any one of the four intersecting grid).

As mentioned, solution in this case may not be unique, what is expected here any of the all possible solution where count of grids in the solution is minimized.


texts: Understanding the problem of applying A* on a multilevel network

texts:
I am currently working on an indoor routing project for a multilevel building.For this, I want to find a shortest path over multiple level (or 'floors'). I have read some papers on A* in such an environment and I think I understand why it is needed to first let the algorithm calculate the shortest path to the next structure to overcome the vertical movement. As an example, consider two level of which the source node is on the first floor and the destination on the second. First the algorithm should calculate the path to the closest elevator in the first floor. However, I only THINK that I have understood it, but I cannot find any proof of it if I'm really correct. Maybe one of you could have a short look at the following image. This image shows the case where I think A* would not find a path as the vertical movement structure (e.g. the elevator, which would be here the grey nodes) is not set as a destination. A* would 'heat in the wrong direction' and ignore the nodes which actually would lead to the destination. Here, in the image, it ignores the outgreyed nodes which are representing the elevator. Is this a correct example for the difficulty you face when applying A* on a multilevel network?

I understand, that the used heuristic (here Euclidean distance) would lead the algorithm to, roughly spoken, 'search in the wrong space' in which it cannot find the destination node. Further, the heuristic prevents to search in other locations, like a barrier.


texts: Can I set individual scale levels for an ArcGISDynamicMapServiceLayer?

texts:
I have a map service that contains 6 layers.

I turned one layer into a feature layer. I used this layer in a query and show the result.

I then would like to just show the other 5 layers as a ArcGISDynamicMapServiceLayer, but I didn't set any min/max scale levels for the layers when I published the MXD. I figured it'd be best to do this programmatically, but I can't figure out how.

Can I set the min/max scale levels independently for each layer in the ArcGISDynamicMapServiceLayer? If not, what's the best way to approach this?

Thanks!

EDIT

Can someone explain to me how this question isn't clear? I don't believe it could get any more clear.


texts: Can I take an intersection for polygons WITH A COMMON ID?

texts:
I have a thiessan polygons for a bunch of cell towers, and need to cut them when they extend more than 12 miles from a cell tower. What I was thinking of doing is taking a buffer from each point, and taking an intersection of each tower's thiessan polygon and its buffer, but I don't know how to take an intersection BY A UNIQUE ID (the tower name). Can anyone help?

(I have ARCinfo 10, for what it's worth)

Thanks!

Nick


texts: Finding closest marker to user location Leaflet KNN

texts:
I am using Leaflet KNN and would like the closest station to a users location, from the web browser, to have a pop up saying it is their closest station.
How can I find the nearest marker to the user location?
Using KNN I am currently able to find the closest station to the coordinates of a users mouse click, I am also able to find the users location but not the closest station to the users location.
//User location
var Position, Accuracy;
function onLocationFound(e) {
  if (Position) {
    Leafmap.removeLayer(Position);
    Leafmap.removeLayer(Accuracy);
  }

  var radius = e.accuracy / 2;
  //pop up stating the users position
  Position = L.marker(e.latlng)
    .addTo(Leafmap)
    .bindPopup(
      "
You are here

You are  within " +
        radius +
        " metres from here
"
    )
    .openPopup();
  Accuracy = L.circle(e.latlng, radius).addTo(Leafmap);
}

//error message if user location cannot be found
function onLocationError(e) {
  alert(e.message);
}
//display the layers to the map
Leafmap.on("locationfound", onLocationFound);
Leafmap.on("locationerror", onLocationError);

//locate function
function locate() {
  Leafmap.locate();
}
// call locate every 10 seconds... forever
setInterval(locate, 10000);

content_index = leafletKnn(contentLayer);

Leafmap.fitBounds(contentLayer.getBounds());

Leafmap.on("click", function (ev) {
  var nearestResult = content_index.nearest(ev.latlng, 1)[0];

  nearestResult.layer.bindPopup("I am your closest station!").openPopup();
});


texts: I'm new to digitization. I have the following map and I'm trying to digitize the purple features. Is there any way to automatically digitize all the purple lines on the map.


texts:
Recently started using QGIS 2.18, for basic mapping.
I am using Swedish GEO-agancy data shapefiles as bottom layers with my own shapefiles as toplayers (EPSG:3006 (SWEREF99 TM)), including a GPX file (EPSG:4326 (WGS84)) manipulated by GPXeditor.

When using the Print Composer to export the map as a GeoTIFF I can not choose a "Reference Map" under "Composition -> Export Settings".

So the exported TIF does not include any georeferencing when opening it in ex. Avenza on my phone.
Also the map on my phone is already 2/3 zoomed in, I can't see an overview by zooming out. I can zoom in though.

All coordinates match up in the QGIS project window with the EPSG:3006 system, even the GPX waypoints that are Google Maps system coordinated EPSG:4326 (WGS84).

Why can't I set a Reference Map ?
Is it because all layers are EPSG3006 exept the GPX waypoints that are EPSG:4326?

Or is there some other reasone the map does not export a GeoTIFF.


texts:
Consider a region (2D) filled with lines randomly (following Figure). We are interested in filling the empty spaces between lines including four boundary edges in a way:

0- maximizing the size of parcels;
1- shape of filling parcels is square aligned horizontally or vertically;
2- shape of filling parcels is square, i.e., relaxed alignment;
3-  shape of filling parcels is any quadrangle. our original question

So for now there are three different scenarios.
Note that the lines are of the form [x1,y1,x2,y2] point set, real numbers.

[* * *] Ideas of possible solutions/algorithms/code snippets/etc are more than welcome.

Update 1: We could manage a solution for the first case:

Steps are:
1- lines
2- rasterizing lines into a bitmap
3- searching nearby cells for each cell of desired color (i.e., the same color) with an objective function to maximize the area i.e., the number of cells.

It works well however it covers only the first scenario and it is also slow.

Update 2:
We assumed that the reader is familiar with the concept of space-filling-tiling. You may follow the link for inspiration. However note that our problem is different. As we do not fill the empty space randomly and we do not chose the size randomly. The solution should be iterative. For all the cases, there is no limit on the number of parcel being fitted. Indeed, it is up to user to limit the iteration number, by choosing a minimum area for parcels, for example. This is obvious in the example given above in which we discretized lines into pixels with specified size. That is, the procedure should run until entire empty area is filled respecting the criterion e.g., the maximum area of parcels.

Update 3:
summary:
One application is to find out the distribution of extractable intact 'rock' blocks in a heavily fractured 'mine'. This could be very helpful for many aspects including drilling design, financial evaluation etc.
description:
For a mine of decorative rock (stone) the products which are the blocks of intact rocks cut as rectangular cubes the price is closely dependent to the size of the block.  Extraction of a block from a suitable area i.e., with no major fracture will be desired if the amount of remaining parts is small as possible. Usually, the small pieces of rocks have no economic value relatively and are considered so as waste.
The question in this post investigates solutions for this kind of problem.

A mathematical view for the problem can be stated as follows:
2D: Find all rectangles that could be extracted from a given 2D region with some lines optimized for larger rectangle size as possible.
3D: Find all rectangular cubes that could be extracted from a given 3D region with some sub-planes (better: polygons) optimized for larger block size as possible.

Since this is part of an on-going research, some of the questions asked in the comments below do not have certain answers that we can provide. We believe that the provided information here so far is indeed enough to get the overall picture of the problem.  Nevertheless, we provide some details as we can for community benefits.
You may put some restrictions on the solution for the ultimate question, although, we believe it is always possible to add more later.  For example, follow these:
{2D case}
The best size of a block (economically optimum rectangle) to be extracted under the conditions mentioned above, is 1x1 m given 10x10 m for the region in the example. This is one constraint defined based on economical value. The minimum workable size for cutting etc, let be 0.15x0.15 m; so this the second size limit.

The figure above shows the economic value function depending to the block size. So for this particular case every rock piece smaller than 0.15x0.15 m is just waste. There will be no block size larger than 1.7x1.7 m due to operation limits.


texts: Obtain CRS from existing axis with cartopy

texts:
Is there a method to get the CRS from a figure created with cartopy?
Say for example, I only had access to the ax created here:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 3))
ax = plt.axes(projection=ccrs.PlateCarree())

I was hoping to be able to do something like ax.get_crs() to return 'PlateCarree'


texts: How to deal with two dimensional lat and lon variables in a NetCDF file in R?

texts:
I'm trying to convert a NetCDF file to a dataframe in R using the RNetCDF package. The NetCDF file is very different from what I have experience with, in that the lat and lon values each have two dimensions:

  data<-open.nc("myfile.nc")

  lon<-var.get.nc(data,"lon")

  dim(lon)

[1] 232 196

and also when I inquire about "lon":

  var.inq.nc(data,"lon")

$id

[1] 1

$name

[1] "lon"

$type

[1] "NC_DOUBLE"

$ndims

[1] 2

$dimids

[1] 1 0

$natts

[1] 4

and here is the summary related to lat and lon:

  print.nc(data)

dimensions:

    south_north = 196 ;

    west_east = 232 ;

    time = UNLIMITED ; // (2920 currently)

variables:

    double lat(west_east, south_north) ;

            lat:units = "degrees_north" ;

            lat:long_name = "latitude coordinate" ;

            lat:standard_name = "latitude" ;

            lat:description = "Latitude, South is negative" ;

    double lon(west_east, south_north) ;

            lon:units = "degrees_east" ;

            lon:long_name = "longitude coordinate" ;

            lon:standard_name = "longitude" ;

            lon:description = "Longitude, West is negative" ;

I don't understand what the two dimensions mean and how I can convert them into a vector and use it to create a data frame. Has anyone had experience with a data set like this? What is a convenient way to deal with it in R?


texts:
I am trying to mask out some of the Sentinel-2 L1C Bands prior to other applications that require me to keep the classic .SAFE structure of the data.
I usually use rasterio with no trouble, but even with an operation as simple as
import rasterio as rio

input_band = ".../T33NVE_20180101T091351_B02.jp2"
out = ".../test"

with rio.open(input_band) as src:
    data = src.read()
    prof = src.profile

with rio.open(out, mode='w', **prof) as dst:
    dst.write(data)

Here is what I get :
Input Image :

Output Image

I know JP2000 is a lossy format, and I don't usually use it....
What can I do to end up with the jp2 files recognized by Sen2cor?


texts:
I'm having a bit of trouble with the Alias variable in MapBasic.
In short, I have a table where I was to group by a specific column.
Then I want to run through each unique value from the grouped table
through the original one and combine each one.
So I have a table of cities with a column of which country each city is in.
My program gives the user to choose a table and a column within that table to gorup by.
Then it basically goes like this:

Dim aSelection, aGrouped as Alias

Selection = selectedTable & "." & selectedColumn

Select * from selectedTable into Temp_GroupedTable Group By selectedColumn

aGrouped = Temp_GroupedTable & "." & selectedColumn

So a) is this code correct? b) How can I fetch the data from each row in the
Temp_GroupedTable so that I can compare it with the values in the original table?
I currently have this general setup within a loop:

Fetch rec i from Temp_GroupedTable

Select * from selectedTable into New_Temp_TAB Where aSelection = aGrouped

Objects Combine

I know there is something wrong... But I don't know how to fix it.
Any ideas?


texts:
I have a model which adds and calculates a new field with a workspace path.
I upload only one layer from database to the model (called "tline"). Parse path can detect only the name and file of this layer, not the whole path to the GDB file which contains this layer. The path can be found only if I add a workspace variable and set the path to database in this variable. So when I set a parameter "Parse type" like "Path" for "tline" layer, model won't work, it just add a new field and says that "The process did not execute because the precondition is false."

How can I find the path to GDB automatically using model builder? It is also important when I have several GDB files that have the same layer "tline" within. I just want to overcome a manual input of GDB file into the model.


texts: Is it possible to do a CQL filtering on an array field in a Postgis DB

texts:
I have a Postgis (2.2.3 in PostgreSQL 9.5) database with a view that has aggregate array fields.

I can filter on these fields directly on the DB with

Select * from myschema."ViewWithArrays" v
where '7415CM' = ANY(v."Postcodes")

But I cannot seem to find the syntax for filtering these array fields with CQL to my Geoserver (2.10.1) which has the view as a layer.
http://docs.geoserver.org/stable/en/user/filter/ecql_reference.html#filter-ecql-reference

The IN keyword seems logical, but the following does not work:

CQL_FILTER=7415CM IN Postcodes

This gives the following error:

  Could not parse CQL filter list. Encountered "7415 CM" at line 1,
  column 25. Was expecting one of: "not" ... "include" ... "exclude" ...
  "(" ... "[" ... "id" ... "in" ...  ...  ...
   ... "-" ...  ...  "*" ...
   "/" ...  "+" ...
  "-" ...  "not" ...  "like" ...
   "ilike" ...  "exists" ...
   "does-not-exist" ...  "is" ...
   "between" ...  "=" ...
   ">" ...  "<" ...
  ">=" ...  "<=" ...  "<>" ... Parsing
  : VerzamelingId = 241 AND 7415CM IN Postcodes.

if you need the stacktrace of the geoserver, let me know, I will add it than

edit:

putting the statement 7415CM in single quotes does not help:

  Could not parse CQL filter list. Encountered "\'7415CM\' IN" at line
  1, column 25. Was expecting one of: "not" ... "include" ... "exclude"
  ... "(" ... "[" ... "id" ... "in" ...  ...  ...
   ... "-" ...  ...  ...
   ...  "*" ...  "/" ...
   "+" ...  "-" ...
  "not" ...  "like" ...  "ilike" ...
   "exists" ...  "does-not-exist" ...
   "is" ...  "between" ...
   "=" ...  ">" ...  "<"
  ...  ">=" ...  "<=" ...
   "<>" ... Parsing : VerzamelingId = 241 AND '7415CM'
  IN Postcodes.


texts:
The below code exports monthly data. However while it is saving in Google Drive, returns this error:
Error: Expected a homogeneous image collection, but an image with incompatible bands was encountered: First image type: 1 bands ([1993030100]). Current image type: 1 bands ([1993030200]). Image ID: null Some bands might require explicit casts.

code link: https://code.earthengine.google.com/a3062bac73d72ca444106081ff99ecca
var temporalCollection = function(collection, start, count, interval, units){

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

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

    var startDate = originDate.advance(ee.Number(interval).multiply(i), units);
    var endDate = originDate.advance(ee.Number(interval).multiply(ee.Number(i).add(1)), units);

    return collection.filterDate(startDate, endDate).mean().clip(table)
    .set('system:time_start', startDate.millis())
    .set('system:time_end', endDate.millis());
  }));
};

var sailinity = function(img){
  var nam = img.id();
  var bands = img.select(['salinity_0']).multiply(0.001).add(20.0).clip(table);
  return bands.rename(nam)
  .copyProperties(img,['system:time_start','system:time_end']);
};

var HYCOM = ee.ImageCollection("HYCOM/sea_temp_salinity")
.filterDate('1993-01-01','2003-01-01')
.filterBounds(table)
.map(sailinity);

print(HYCOM);

var HYCOM9303 = temporalCollection(HYCOM, '1993-01-01', 120, 1, 'month').toBands();

print('HYCOM9303',HYCOM9303);

Export.image.toDrive({
  image: HYCOM9303,
  description: 'HYCOM9303',
  scale: 10000,
  region: table,
  crs: HYCOM9303.getInfo().crs,
  maxPixels: 1e13
});


texts: How to add more than one xml file per layer

texts:
I'm using this code to add xml (osm relations) files to my map.

osm_layer = new OpenLayers.Layer.GML("Polygon", relation.xml, {
format: OpenLayers.Format.OSM,
style: {strokeColor: "blue"},
projection: new OpenLayers.Projection("EPSG:4326")
});

But what I need is to add more than one relation.xml per layer, as relation1.xml, relation2.xml in the SAME layer. Until know I can only show multiple xml if they go in multiple layers.

When using vectors I just created the layer and then added the vectors by something like

polygonFeature = new OpenLayers.Feature.Vector(
                 new OpenLayers.Geometry.Polygon(linearRing), style);

And then

Layer.addFeatures(polygonFeature);

However, I can't make it work with xml / OSM relation format...

Thanks!


texts: StyleFunction Openlayers-3 - multiple callback variables?

texts:
I am using the "styleFunction" from the following example: http://openlayers.org/en/v3.9.0/examples/geojson.html

var styleFunction = function(feature, resolution){
    var properties = feature.getProperties();
    var type = feature.getGeometry().getType();
    switch (type) {
        case 'Point':
        ...return some styling...
        break;
        case 'Polygon':
        ...return some styling...
        break;
    }
}

I am currently using only the "feature" from the callback. I can't seem to find documentation about this function on Openlayers.org so I ask here:

Are there more variables available for callback and is there a way to find out without documentation?

(I need to know the layer name the function is styling, and I can't get that from the feature alone.)


texts: Intersections between two Linestrings - SpatiaLite

texts:
I have two LineStrings and would like to draw the point of their intersection.
The two lines are similar to this:

LINESTRING(445798.43 6153379.87, 445892.89 6153340.08, 445987.35 6153300.3, 446081.82 6153260.51,
446176.28 6153220.73, 446270.75 6153180.95, 446365.21 6153141.16, 446459.67 6153101.38, 446554.14
6153061.59, 446648.6 6153021.81, 446743.07 6152982.02, 446837.53 6152942.24, 446931.99 6152902.46,
447026.46 6152862.67, 447120.92 6152822.89)'

I am sure that the two lines are not parallel, but when I execute this query:

SELECT ST_Intersection(p.geometry, r.geometry) AS geom FROM tableOne AS p ,
table2 as r WHERE st_intersects(p.geometry, r.geometry)

I get a return value equal to NULL.

I don't know if it can be relevant, but if I run the following query:

select geometry from tableOne

It return: BLOB sz=11104 GEOMETRY

Can anyone help me? thanks


texts: Openlayers 3 + WFS GML2 layer wrong projection and displaying on wrong location

texts:
I have requested vector layer through WFS GetFeature from GeoServer in GML2 format but when displaying WFS layer its showing on different place. same layer i have requested in format Geojson its showing on correct location. My vector layer is in EPSG:4326 projection and when displaying in Geojson format its correctly displaying but in GML2 format its showing on wrong place which is on EPSG:900913 location. As i understand WFS layer is in EPSG:4326 projection in Geoserver but whenever layer requested in GML2 format its converting or changing projection to EPSG 900913 and showing on wrong place. My question is how to display GML2 layer in correct position(EPSG4326) and why GML2 layer automatically changing projection? or if GML2 showing on EPSG900913 what is the code i can apply so it will be in EPSG4326 ?
I am giving here my code also:https://jsfiddle.net/seema30/dww9gxcy/ please give the suggestions.


texts: osm2pgrouting: Import multiple OSM files?

texts:
I want to import multiple OSM data files for a few different countries into postgres using osm2pgrouting.

Can I just run osm2pgrouting on each file with the same database settings (e.g., import countryA.osm, then countryB.osm, etc.) or would this cause problems?


texts:
I have recently installed the Google Maps Engine Plugin for QGIS 2.0 and it doesn't seem to work. Whenever I click on the Log in button it pops up with:

Google
Error:invalid_client
Bad request.

cookie_policy_enforce=false
response_type=code
scope=https://www.googleapis.com/auth/mapsengine https://www.googleapis.com/auth/userinfo.profile
redirect_uri=urn:ietf:wg:oauth:2.0:oob
approval_prompt=auto
state=123456789
client_id=myemailaddress@googlemail.com

I have uninstalled and re-installed the plugin several times but it doesn't seem to make a difference. Any ideas would be appreciated.


texts: JS libraries that work with land use maps?

texts:
I'm incredibly new to GIS, but not to web development.

I have two land cover maps for two islands, which I'm spatializing data on top of for a web app. This data is agnostic of any real location data - I know that certain data points belong on certain land types, but since everything is just images and arrays, I've never needed to worry about long/lat. Up until now I've been using p5.js to build out the maps and the points.

Someone accidentally mentioning GIS triggered a quest to see if there are any JS libraries that deal with data from land use maps. Whether they generate GeoJSON from land use/cover images, or can somehow intake land cover data and display it over top of a live map, I'd love to see if I can future-proof my web app by using real maps instead of screenshots off of Google Earth. "Land cover/use" doesn't seem like it's a frequent tag on these libraries, but I'm new to this world and might be missing something.

Any ideas?


texts: Provide MbTiles from GeoServer with mutliple vector Layers

texts:
Is there a way to serve pbfs tiles by loading a mbtiles file (multiple vector layers) into the GeoServer?

For this we used before Klokantech TileServer Gl to offer big mbtiles (osm vector tiles) to load it into mapbox-gl.
Now we want to move to GeoServer, and the question is if we may only use GeoServer.

I find different GeoServer extentions and there are some question marks:
- What is the GeoServer MBTiles Data Store exactly? only Images, or also vector data?
- Does Vector Tiles Extension also support multiple layers? And will it also work with mbtiles file?


texts:
I try to do something simple, use official data of Belgium zones that I found here: https://eservices.minfin.fgov.be/myminfin-web/pages/cadastral-plans.

So I downloaded the files of one city and opened it in QGIS. That works but the coordinate system seems to be totally wrong. I would like to switch to an international standard.

In fact, I need to use the Googlemap API to get coordinates from full addresses and use that to get the position. Is there a way to correct that?


texts: Losing resolution using gdalwarp to change projection?

texts:
I’m trying to project an image from Mars using gdalwarp, from eqc to polar stereo projection:

gdalwarp -overwrite
-s_srs "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs"
-t_srs "+proj=stere +lat_0=-90 +lat_ts=-90 +lon_0=0 +k=0.1 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs"
-of GTiff MarsEqc.TIF MarsStere.TIF

However, in the result I lost spatial resolution, as you can see in the image. As the original resolution is 1km/px I tried to change -tr and -ts options but with no good results.

If I change the projection directly in QGIS the result is pretty good, without losing resolution. Of course to use QGIS would be a solution, but I have to process a lot of images.

What parameters do I have to add to maintain the original resolution?


texts: How to create a WMS Service to display a geo-reference image

texts: I currently have little experience in the opensource area and I'm guessing this is where I need to go if I want to create a WMS Service without needing ArcGIS Server.  Basically my goal is to take a local geo-referenced image (geo-jpeg) and create a WMS service to display that image.  I've done some searching on this topic and I can't seem to find any specific information on this.  I would prefer not needing to create a MapServer installation but if that's one of my best options, I will try going that route.  So at this point I was hoping for some direction on how to create a WMS service, hopefully without needing to setup my own server.  Maybe a blog post or tutorial on doing this.  Thanks...


texts: raster2pgsql missing in Opengeo-suite installation

texts:
I installed opengeo-suite 3.0.2 on Ubuntu 12.04 and although the stack seems to be working fine, raster2pgsql is no where to be seen!

root@mywebsite:/usr/local/lib/postgresql/9.2# locate raster2pgsql
/usr/lib/postgresql/9.2/bin/raster2pgsql

root@mywebsite:/usr/lib/postgresql/9.2/bin/raster2pgsql# raster2pgsql
/usr/bin/raster2pgsql: No such file or directory

Did the installation go wrong? Or is it really not included in the opengeo-suite?

Update

ls -l /usr/bin/raster2pgsql
ls: cannot access /usr/bin/raster2pgsql: No such file or directory

ls -l /usr/lib/postgresql/9.2/bin/raster2pgsql
-rwxr-xr-x 1 root root 159016 Sep 24  2012 /usr/lib/postgresql/9.2/bin/raster2pgsql


texts: Calculate traveltime from a Least-Cost Path route

texts:
I'm working in least-cost paths. I used r.walk to create cost maps. I used SAGA least-cost path routine to create routes from several points to the cost-map center.

Now I need to calculate the travel time of these routes. How I can do it?


texts: Check if a land parcel (polygon) has access on a road or not

texts:
I have a feature layer with land parcels as polygons. Polygons (land parcels) are either adjacent in all its around sides by other polygons but some of them have a road a which is not actually a feature layer representing roads but is just an empty space between parcels. In this case, parcels have access to roads. The question is how can I automattically check if a polygon has access onto a road using VBA and AcrObjects.

Thanks
Deme


texts: GDAL not accepting *.tif to list through files in directory

texts:
I am using gdal_merge.py to mosaic tif. If I explicitly write the file names it works, but if I use *.tif it fails with the following error.

Traceback (most recent call last):
  File "C:\Program Files\GDAL\gdal_merge.py", line 509, in <module>
    sys.exit(main())
  File "C:\Program Files\GDAL\gdal_merge.py", line 392, in main
    ulx = file_infos[0].ulx
IndexError: list index out of range

Can anyone tell me why?


texts: GIS Cloud Publisher for ArcMap extension does not work

texts:
GIS Cloud Publisher is an extension from GISCLOUD.COM that allows one to publish from ArcMap to GIS cloud directly. I followed the instruction provided on the website concerning the installation but it does not work. I does not appear in the Toolbars in Arcmap.

I am using ArcGIS 10.3

Can anyone help with how to make the the GIS Cloud Publisher extension appear or work in ARCMap.


texts: KEA File Conversion using gdal_translate

texts:
I was using ARCSI (Atmospheric and Radiometric Correction of Satellite Imagery) to atmospherically correct data, the default output format is KEA ('kea'). I'm trying to convert these to 'tif' format using this gdal statement (gdal_translate -of GTIFF LS8_sref.kea reflectance.tif), but was getting an error message (too many command options).

Are there alternative ways to convert kea to tif or any raster format?

@AndreJ - This is the file:https://www.dropbox.com/s/kdy80qggdgufum5/LS8_20140927_lat50968500lon5287500_r28p19_rad_sref.kea?dl=0


texts:
I'm reprojecting GFS GRIB files to EPSG:3857 using gdalwarp. Resulting file get little bit shorter to the right so when I zoom in from zoom level 5 and higher, I can see vertical line with missing data. Check image below:

When I check edges with gdalinfo I could clearly see that map is indeed shorter:

Upper Left  (-20037508.343,26491625.382)  (180d 0' 0.00"W, 88d12' 0.00"N)
Lower Left  (-20037508.343,-26492654.355) (180d 0' 0.00"W, 88d12' 1.05"S)
Upper Right (20036608.125,26491625.382)   (179d59'30.89"E, 88d12' 0.00"N)
Lower Right (20036608.125,-26492654.355)  (179d59'30.89"E, 88d12' 1.05"S)
Center      (    -450.109,    -514.486)   (0d 0'14.56"W,  0d 0'16.64"S)

For tile x = 32767; y = 14181; zoomLevel = 15; I calculated geo coords:

x0 = 20036285.349764004
x1 = 20037508.34221653
y0 = 2693029.380466357
y1 = 2694252.372918885

and I get complexly blank image.

For tile x = 63; y = 28; zoomLevel = 6; There is missing data vertical line at the right side of image.

I am using SOURCE_EXTRA=1000. I tried much higher numbers, but I guess resulting tif bounding box is not calculated correctly so this doesn't help.

Downloaded any grib2 file from here

Extracted 270 variable using wgrib2 tool with following grid params:

-new_grid latlon 180:23040:0.015625 -90:11520:0.015625

Finally, I used gdalwarp command:
gdalwarp -t_srs EPSG:3857 -wo SOURCE_EXTRA=1000 --config CENTER_LONG 0 -of GTiff --config GRIB_NORMALIZE_UNITS NO -co "TILED=YES" 270-16.grb test.tif

Extracted and interpolated grib file: link

Resulting TIF file: link

Any thoughts how to fix this?


texts: QGIS- Widget Photo Attachment do not show photo (small pix)

texts:
I have attached photos with their file path in the column. And I have set the height and width of the photo. but it does not work to show the photo on the attribute table.

and this is when I show the attribute table with Identify Features icon in QGIS. but the photo that I attached on attribute table just show in a small size although I set photo height and width wit 500px each. you can see the screenshot below.

Can any one fix this problem?


texts: Is it possible for a script to be built that will use a reference point and have it draw several polygons based on distance and bearing, then assess the terrain under the polygon for highest obstacles or terrain?  I know there is a gent on github that has the formulas built to identify the locations, I need to extend that capability by assessing underlying terrain.


texts:
Is it possible to adjust tabel headers of GeoJSON files in ArcGIS?

If yes, which max length has a row header?

In QGIS there is no possibility to edit a GeoJSON table.


texts:
In a previous question (convert-xyz-points-to-topography-and-calculate-mean-slope), I asked how to interpolate point elevation data to a raster in order to measure slope. Thus I use the Raster > Interpolation tool, and create the raster with TIN method, grid 300x300.

This works, however when I calculate slope (raster->Analysis -> DEM) it gets awful.

I have few points, but I know some softwares such as Surfer provide tools to obtain a quite smooth surface by resampling.

Is that possible in QGIS? Do I need to resample the created raster? Or is there another way to get more "continuous" slope values ?


texts:
How to switch on snap option in QGIS 2.18 when creating customize polygon as in:

class customPoly(QgsMapToolEmitPoint):

def __init__(self, iface):
    self.iface = iface
    QgsMapToolEmitPoint.__init__(self, self.iface.mapCanvas())
    self.rubber=None
    self.dup_Pt=None
    self.cnt=0
    .....

def canvasPressEvent(self, event):
    if event.button()==Qt.LeftButton:
        if self.draggable:
            self.pos=event.pos()
            mapping_Pt, lyrPt=self.transformCoordinates(self.pos)
            ....
            self.rubber.addPoint(mapping_Pt,self.cnt+1)
            self.dup_Pt.addPoint(mapPt,self.cnt+1)

            ....

The above code is part of custom polygon, I want to know, how to switch on snapping while picking point or QgsMapToolEmitPoint or in canvasPressEvent(self, event) or any other best option to switch on snap while in pick point while creating custom polygon


texts: Display markers and autohide them after few seconds using data of nested JSON object from an API in Leaflet Map

texts:
I'm trying to load 4 markers at a time and timeout should be 5 secs for them and then new markers should come up and old one's should be hidden automatically.
I have written this code and unable to proceed further. I have tried out Setinterval() and SetTimeout() but both of these don't seem to work.
I actualy want to display markers and hide them on the basis of "date" which is a nested JSON Object but I don't know how to.

JSON data from an API:

   {"datapoint":[{"id":"H0031","In_time":{"date":"2017-10-24 9:41:17.783000"},"lat1":"24.55","long1":"71.75","lat2":"27.715186","long2":"83.323999"}]}
{"datapoint":[{"id":"H0032","In_time":{"date":"2017-10-25 10:41:17.783000"},"lat1":"24.55","long1":"71.75","lat2":"27.715186","long2":"83.323999"}]}
{"datapoint":[{"id":"H0033","In_time":{"date":"2017-10-26 11:41:17.783000"},"lat1":"24.55","long1":"71.75","lat2":"20.715186","long2":"89.323999"}]}
{"datapoint":[{"id":"H0034","In_time":{"date":"2017-10-27 12:41:17.783000"},"lat1":"21.55","long1":"61.75","lat2":"29.715186","long2":"81.323999"}]}
{"datapoint":[{"id":"H0035","In_time":{"date":"2017-10-28 01:41:17.783000"},"lat1":"28.55","long1":"77.75","lat2":"27.715186","long2":"83.323999"}]}

Code:

fetch('http://www.example.com/?qty=10')
        .then((response) => response.json())
        .then((response) => {
            let datapoint = response.datapoint;
            for (let i=0; i<datapoint.length; i++) {

                let hlat = parseFloat(datapoint[i]["lat1"]);
                let hlon = parseFloat(datapoint[i]["long1"]);
                let glat = parseFloat(datapoint[i]["lat2"]);
                let glon = parseFloat(datapoint[i]["long2"]);
                let markerLocation1 = new L.LatLng(hlat, hlon);
                let marker1 = new L.marker(markerLocation1,{icon: greenIcon});
                marker1.addTo(map).bindPopup(popup,customOptions);
                let markerLocation2 = new L.LatLng(glat, glon);
                let marker2 = new L.marker(markerLocation2,{icon: BlueIcon});
                marker2.addTo(map).bindPopup(popup,customOptions);
                var latlngs=[];
                let latlng1=[hlat,hlon];
                let latlng2=[glat,glon];
                var offsetX=latlng2[1]-latlng1[1];
                var offsetY=latlng2[0]-latlng1[0];
                var r=Math.sqrt(Math.pow(offsetX,2)+Math.pow(offsetY,2)), theta = Math.atan2(offsetY,offsetX);
                var thetaoffset=(3.14/10);
                var r2 = (r/2)/(Math.cos(thetaoffset)),
                    theta2 = theta + thetaoffset;
                var midpointX = (r2 * Math.cos(theta2)) + latlng1[1],
                    midpointY = (r2 * Math.sin(theta2)) + latlng1[0];

                var midpointLatLng = [midpointY, midpointX];

                latlngs.push(latlng1, midpointLatLng, latlng2);
                var pathOne = L.curve(['M',latlng2,'Q',midpointLatLng,latlng1], {snakeArray: 5, color:'grey' ,animate: {duration: 2500, iterations: 2}}).addTo(map);

            }

        })

P.S- I have connected one marker to the corresponding one using L.curve function which generates a curved line.


texts: Adding field vector layer in QGIS modeler and using it in next steps?

texts:
I am trying to create a model with QGIS to intersect a polygon layer P with a vector grid G and find the % of P for every cells of G. These are the steps to find this :

create grid G with id unique field
find GRIDSURF = cell area
intersect grid G with polygon layer P
dissolve intersected layer by attribute grid id
calculate PSURF = dissolved polygons area / GRIDSURF

My problem with the modeler is that, in step 5, the modeler can't find the field GRIDSURF created in step 2

What I really don't understand, is that in the following model with 2 outputs, the GRIDSURF field (with $area formula) is in the result_grid output layer, but this same field is not in the layer submitted to INTERSECT ... and I have an error after the DISSOLVE saying that the GRIDSURF field is not found.

Did I miss something ?


texts: How to assign (or project) CRS to LiDAR file in lidR?

texts:
Can someone enlighten me if there is a difference between initializing a CRS object using the EPSG code versus just the +proj string read from a shapefile?

Consider the following objects:

A SpatialPolygonsDataFrame read in from a .shp using readOGR, no additional arguments. I was handed this data.

> itcs[[1]]@proj4string
CRS arguments:
 +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

A lidR tile:

tile<-readLAS(paste("../data/2015/cropped_",fname,sep=""))
tile@crs<-CRS("+init=epsg:32617")
> tile@crs
CRS arguments:
 +init=epsg:32617 +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

To me, these should have the same coordinate system, that is, they have the same projection, zone, datum, ellipsoid and units. One already had a coordinate system, one was assigned using the EPSG code. When plotted they appear to be misprojected, overlapping, but skewed and non-aligned. I'm trying to track down the source of this misprojection.

Perhaps similar here why is there a difference between using EPSG and using basic code


texts: Exact same code gives different result in server and local host

texts:
I have a back end that uses PostGIS to filter points inside a polygon. Backend is in Django, and here is the code that does the filtering:

def get_queryset(self):
    polystr = self.request.query_params.get('poly', None)
    if polystr:
        poly = GEOSGeometry('SRID=4326;' + polystr)
        qs = PropertyPost.objects.filter(location__contained=poly)
        return qs

So using that I would be able to filter the points inside a polygon, like this in local host (Windows):

http://127.0.0.1:8000/prop-posts/?poly=POLYGON((50.844738 28.957921,50.844906 28.959785,50.8433 28.957801,50.844738 28.957921))

the above returns the points inside the polygon.
Now I have 3 points inside the polygon that correctly will be returned by above request.
However I have uploaded my code in digital ocean server on Linux, and the below code returns only two points:

http://127.137.139.1/prop-posts/?poly=POLYGON((50.844738 28.957921,50.844906 28.959785,50.8433 28.957801,50.844738 28.957921))

I am puzzled that what can cause this. The code
is exactly the same. I usually test things locally and push it to server and everything has been pretty similar so far. So it is odd that it is not working.

Both PostgreSQL (Windows and Linux) are 10.6. However PostGIS versions are a bit different (but I doubt this should make such an issue):

POSTGIS="2.4.3 r16312" on Linux server

POSTGIS="2.4.4 r16526" on local host Windows.

Is there some setting that might be different on Windows and Linux?


texts: Combining linesegments that are contionous into one feature in QGIS 3.10.7

texts:
I have a dissolved line, where I have clipped some segments where I would like to create a feature for each individual segment.
I have tried using "Vector Geometry -> Multipart to singlepart", but the linesegments gets divided into even smaller parts, as shown in the figure below.
Is there a way to combine the linesegments into individual features?


texts: Where does "Add XY Coordinates" tool get data from?

texts:
I am using the ArcGIS Add XY Coordinates tool to populate XYZ values in a shapefile after converting .las data to .shp using LAS to Multipoint tool.

Does the Add XY Coordinates tool get the location info from the LAS dataset OR is it getting it from the data frame somehow?

I need to make sure I am getting the actual data from the LAS files because, I will eventually convert the shapefiles dbf to a .pts file to do some other calcs outside of ArcGIS.

Also, is there a better way of getting XYZ from LAS in ArcGIS?


texts: Best projection for Getis Ord Gi* Statistic Analysis

texts:
I am about to calculate the Getis Ord Gi* statistic for some rasters whose extents are global.

In choosing a projection, do I want to preserve distance (e.g. equidistant) to obtain accurate results?
What do people recommend as a decent global projection for such an analysis?


texts:
Possible Duplicate:
  How to change Leaflet Map panes layering order (z-index)?

In the Leaflet Layers Control L.Control.Layers() described here, is there a way to control the z-index (layer order) of the overlays?

I'm talking about the overlay portion they show like this:

var overlayMaps = {
"Motorways": motorways,
"Cities": citiesLayer
};

In their example, they have a tile layer and a maker layer. Of course, the marker layer is always on top regardless of which toggle button was clicked when.

In my case, I have two tile layers. The problem is that whichever toggle button I turn on last, THAT layer appears on top.

If both are off and I turn on Cities, then turn on Motorways, Motorways is on top of Cities and partially blocks them out. Is there a way to control this so Cities always appears on top even if Motorways is turned on after it?

I've tried reordering the code to no avail and am aiming to have maybe five overlays. The code is below. Thanks!

var mapAttribution = 'Attribiution',
            mapOptions = {maxZoom: 8, attribution: mapAttribution},
            mexicoMapUrl = 'http://{s}.tiles.mapbox.com/v3/billted98.map-50h5h3pu/{z}/{x}/{y}.png';
            blueMarbleUrl = 'http://{s}.tiles.mapbox.com/v3/mapbox.blue-marble-topo-jul/{z}/{x}/{y}.png';
            blueMarbleGrayUrl = 'http://{s}.tiles.mapbox.com/v3/mapbox.blue-marble-topo-jul-bw/{z}/{x}/{y}.png';
            roadsUrl = 'http://{s}.tiles.mapbox.com/v3/billted98.roads-major-mexico-select/{z}/{x}/{y}.png';
            citiesUrl = 'http://{s}.tiles.mapbox.com/v3/billted98.msm-cities-and-names/{z}/{x}/{y}.png';

    var mex = new L.TileLayer(mexicoMapUrl, mapOptions),
        bm = new L.TileLayer(blueMarbleUrl, mapOptions),
        bmg = new L.TileLayer(blueMarbleGrayUrl, mapOptions),
        roads = new L.TileLayer(roadsUrl, mapOptions),
        cities = new L.TileLayer(citiesUrl, mapOptions);

    // ADD MAP
    var map = new L.Map('map', {
        subdomains: ["a", "b", "c", "d"],
        center: new L.LatLng(23.342256,-102.019043),
        zoom: 5,
        minZoom: 5,
        maxZoom: 8,
        layers: [mex],
    });

    // LAYER CONTROLS
    var baseMaps = {
        "Standard": mex,
        "Terrain": bm,
        "Terrain-Gray": bmg
    };

    var overlayMaps = {
       "Cities": cities,
       "Roads": roads
    };

    layersControl = new L.Control.Layers(baseMaps, overlayMaps);

    map.addControl(layersControl);


texts: QGIS Raster Calculator dialog doesn't resemble that in documentation

texts: I'm very new to QGIS.  Raster Calculator looks (from online sources) like it ought to have all kinds of buttons and fields.  Mine is rather stripped down.  I have two rasters open, and I just want to divide their values.  Is there a plugin missing or something? (ps -- it doesn't change if I select the rasters)


texts: Get statistics from TIFF and shapefile in QGIS

texts: I have imported to QGIS two files, one .tiff that contains the coverage data in a given area and a shapefile that contains the division of the area (e.g., states). I want to get statistics using both, like for example, I want to know the percentage that each subarea is covered.


texts: ArcGIS API 4.11 JavaScript - "Click" "drag" stopPropagation()

texts:
Using JavaScript in the ArcGIS API, I am querying a feature layer from a Map View, calculating points within a certain radius when the user "click" "drag" with the mouse. I want to be able to stop this functionality when the user presses the spacebar. However, I cannot get this to work. In the JavaScript I turn off event propagation for the "click" "drag" event and then cannot turn it back on again. I understand there is not a "on/off" switch for stopPropagation(). Can anyone recommend a different solution or workaround. An extract of the code I am struggling with below:

var bufferEnabled = false;

function keyDownListener(event) {
  var keyInput = event.key;
  bufferEnabled =
    (keyInput === " " || keyInput === "Spacebar") && !bufferEnabled;
  }

function analysis() {
view.when().then(function() {

  const layer = webmap.layers.getItemAt(0);

  view.whenLayerView(layer).then(function(layerView) {
    watchUtils.whenFalseOnce(layerView, "updating", function(val) {
      // Query layer view statistics as the user clicks
      // or drags the pointer across the view.
      view.on(["click", "drag"], function(event) {
        // disables navigation by pointer drag
        event.stopPropagation();
        queryStatsOnDrag(layerView, event).then(updateCharts);
      });
    });
  });
});

    view.on("key-down", keyDownListener);
    view.on("key-down", analysis);


texts: how to open several MXD files at the same time ( within the same window )in Arcgis10

texts:
I have  a simple question. How can I open several MXD files at the same time in Arcgis 10?
Whenever I try to open the second MXD file Arcgis asks me to save the first MXD file and then it closes it and opens the second one. How can I have these two MXD files open at the same time?


texts: Reading geojson into Leaflet results in error?

texts:
I'm trying to read a geojson file into leaflet. Everything works out when I copy the text from the geojson file and paste it into

var data = { }

So I think the file should be ok for leaflet. But I rather want to be able to read the actual file. Right now I keep all files locally on my hard drive. I loaded Ajax and tried to read the file by

var data = L.GeoJSON.AJAX("data.json")

If I open the HTML the map loads, but not the point markers of the geojson file. When I use the Firefox Console I get:

TypeError: this.callInitHooks is not a function

What am I doing wrong?


texts:
I'm new to QGIS so please bear with me.

I'm using QGIS 3.6.
So I imported a CSV-file with latitude longitude coordinates (WGS 84), created markers and set them to have different colors (to show passing of time) and different sizes (to show precision). The accuracy is given in a column in the CSV file in meters, which need to be represented accurately on the map.

I have set it under Layer properties > Symbolizing > Categorized > Symbol > Size > Fieldtype > diameter (my own variable to show the diameter of the circle I want to draw). The problem is that the diameter isn't accurate. The sizes seem more or less correct in respect to eachother but when I use the measuring tool, the lengths are off. I have tried both to set the units as meters to scale and to set them as map units but this doesn't seem to influence the size. I have also tried to use the Assistant function and selecting the same field but even with different scaling methods I cannot get the diameters right.

So what am I doing wrong? Can anyone help me?