texts:
My project uses a public WMS basemap to fill in data outside my ROI, with the WMS CRS chosen to be suitable given my project CRS.

Now I'm creating an (alternate) layout in my project which outputs in a different CRS (epsg:3857). Unsurprisingly, in this layout the WMS layer has visual artefacts due to reprojection, especially in labels.

Is there a way to change the WMS CRS for the layer to match the layout CRS in that layout only?

As far as I can tell, when I added the WMS layer, my choice of CRS then got hard-coded. I can manually change it by hacking the raster source URI, for instance to crs=EPSG:3857&dpiMode=7&format=image/png&layers=Sub_regional&styles&url=http://maps.geogratis.gc.ca/wms/CBMT (I changed the original CRS= to 3857) and then also setting the layer CRS in QGIS to match.

However, then it looks bad in my canvas in the project CRS, as the problem gets reversed. The WMS offers data in many CRSes, but no one individual one looks good (especially labels) when reprojected to both the project/canvas and layout CRS.

Before I go down some route of having 2 copies of the layer differing only in CRS and URI, one used in this one layout and one everywhere else -- and the resultant hassle of keeping their other layer settings in sync going forward -- is there a less klugey solution? One where QGIS would request the WMS data with a different CRS setting based on the context?

In particular, since a layout can override so many map canvas/project rendering settings, can it modify the WMS CRS, or the layer URI specification more broadly?


texts: Does a large amount of GPU memory help when working with very large images and image mosaics?

texts:
It seems like it would help a great deal to store as much of the image as possible on the GPU for panning and zooming, etc.

We deal with a lot of very large images. Many range from 8 GB to hundreds of GBs in size. At home I've got a 40" 4k monitor and a GTX Titan with 12 GB of video memory, and it performs amazingly well for everything I throw at it. It also looks great when viewing imagery. Is that setup overkill for GIS image processing work?

Would QGIS, ArcMap and PCI Geomatica be able to take advantage of abundant graphics memory when dealing with multi gigabyte images and image mosaics in the range of hundreds of gigabytes?

What GIS applications would see the most benefit from an abundance of GPU memory when dealing with large images?


texts: After changing one layer's projection, layers don't overlay correctly

texts:
I have 2 raster layer in global mapper.
1 layer in WGS 84 and geographic coordinate system.(ALOSDEM)
The second one is in NAD83 UTM Zone.(SASPlanet Image)
Both of these layer are overlayed in global mapper because,as we knew global mapper have the function called "On the fly Projection".

But after i manually change Second layer Using Configure options the result is both of these layer is no longer overlayed.

Why did this happened ??


texts: Displaying density map in ArcGIS Desktop?

texts:
I am using a shapefile from the European Commission on world maritime traffic. They have developed a map of maritime traffic as shown in the comments.

I would like to replicate the same map that I provided with a link in the comments (I am not allowed to use more than 2 links in my post) on ArcMap 10.3 so then I can zoom into Europe and use my map as a discussion for a research paper on the European maritime industry. I have downloaded the point shapefile of the data from the European Commission site. When I add it to my map document, there are lots of uniformly spaced dots as shown here.

So then, I tried to use the Spatial Analyst tool to better depict the densities. I went to Spatial Analysis -> Point Density -> Used the point feature, Output Cell Size: 0.716, Neighborhood: Circle, Neighborhood Settings: Radius: 6, Units: Map. I then made an edit in the Properities -> Symbology tab -> Natural Breaks. Now my map looks like this, but it does not quite look like the map that the European Commission made.

How can I fix this?


texts: GDAL/ogr2ogr ( to convert osm to shapefile)

texts:
I use GDAL / ogr2ogr to convert osm to shapefile but I get an error.

C:\Users\alo>ogr2ogr -f "ESRI Shapefile" outputdir "C:\Users\alo\Desktop\London.osm"

Warning 1: Cannot find osmconf.ini configuration file

ERROR 1: Could not parse configuration file for OSM import
FAILURE:
Unable to open datasource C:\Users\alo\Desktop\London.osm' with the following drivers.
  ->netCDF'
  -> PCIDSK'
  ->JP2OpenJPEG'
  -> PDF'
  ->DB2ODBC'
  -> `ESRI Shapefile'


texts: How to use IGxDialog to save output feature class

texts:
I found this code to take input through IGxDialog box from a guidebook:

IGxDialog gxd = new GxDialogClass();
gxd.AllowMultiSelect = false;
gxd.ButtonCaption = "Add FileGDB";
gxd.Title = "Add All FeatureClasses inside FileGDB";
gxd.RememberLocation = true;
IGxObjectFilter gxObjFilter = new GxFilterFileGeodatabasesClass();
gxd.ObjectFilter = gxObjFilter;
IEnumGxObject gxEnumObj;
gxd.DoModalOpen(ArcMap.Application.hWnd, out gxEnumObj);
IGxObject gxObj = gxEnumObj.Next();
//getting the address of fileGDB
string fileGDBAddress = gxObj.FullName;

I could not figure how save a feature class through IGxDialog Box, there are no examples on internet with arcobject 10.0


texts: GEE - Sentinel 1 - Instrument Mode 'SM' doesn't return any results

texts:
I am trying to visualize some Sentinel 1 strip maps from Tristan da Cunha in GEE (which I know exist based on scihub) - but my code doesnt return them. Any suggestions?

<  // Filter the collection - VV only depending on geometry and date
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('instrumentMode', 'SM'))
    .filterBounds(geometry)
    .filterDate('2019-01-01','2019-07-31')
    .select(['VV']);
print(s1);

//visualisation parameters!
 var viz = {
 bands: ['VV'],
  min: -40.0,
  max: -10,
};

// display each image in collection
function addImage(s1) {
  var id = s1.id;
  var scene = ee.Image(s1.id);
  Map.addLayer(scene, viz, id);
  Map.centerObject(geometry,6);
}

// use map on client-side
s1.evaluate(function(s1) {
  s1.features.map(addImage);
});>


texts:
I'm really interested in knowing what kind of algorithm does GeoServer use when presenting its maps.

GeoServer has its own algorithms, for instant the labelling algorithm. Is there a possibility to change it or modify the algorithm or add something to it? If so where can we find it? How we can change it if we want to?

I see that the algorithm is missing a few standards and it would be helpful if it could be modified.


texts:
These days i am trying to develop web map using FOSS software. I have a WMS in GeoServer but it does not view in ArcGIS Viewer for Flex.

Can you give me XML code for it?

This is my WMS.

http://localhost:8080/geoserver/sac/wms?service=WMS&version=1.1.0&request=GetMap&layers=sac:Luse&styles=&bbox=10.8,1.81,10.57,19.1443&width=512&height=443&srs=EPSG:5235&format=application/openlayers

And is it required crossdomain xml? If yes, how to do it?


texts:
Very new to Leaflet, but learning fast. In the following code example, where/how would I define maxHeight to get the scrollbar in popup for data?

var map = L.map('map', {
    center: [40.75, -74.15],
    zoom: 13
});

var osmAttrib = 'Map data © OpenStreetMap contributors';

var osm = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution: osmAttrib
}).addTo(map);

var cities = new L.LayerGroup([
    L.marker([40.72801, -74.07772]).bindPopup('Jersey City')
]);
cities.addTo(map);

L.imageOverlay(
    'http://research.ccountync.com/maps/images/yellow.png',
    [[40.712216, -74.22655], [40.773941, -74.12544]], {
    opacity: 1
}).addTo(map).bindPopup("This is a yellow box");


texts: May I use OSM in pgRouting to determine TSP solution?

texts:
I have a TSP problem that I'd to determine if my solution path is correct. As stated, a classic TSP, I have a number of locations (bikeshare stations) within a city centre (Glasgow, Scotland) that I'd like to find the lowest cost route in terms of distance which connects all points and arrives back at the start point. There are two sets of locations 22 points and 21 points each.

I've fairly comfortable with QGIS, and would like to take my knowledge a further. I'm considering using pgRouting to uncover the route solution listed above. Right now, I'm confident that pgRoute will output a route. However, I want to confirm that pgRouing can use OpenStreetMap as the network parameters. That the routes produced by pgRouting will follow the road system written in OSM. Specifically, - for example, turns not possible, one-way streets.


texts: Size of output geotiff file not as in the .vrt file

texts:
I have the following vrt file (input.vrt) with GCPs in EPSG:3857

  <VRTDataset rasterXSize="5409" rasterYSize="5408">
  <Metadata />
  <Metadata domain="IMAGE_STRUCTURE">
  <MDI key="INTERLEAVE">BAND</MDI>
  </Metadata>
  <GCPList>
  <GCP Id="" Pixel="678.0000" Line="1154.0000" X="1.632931421600E+07"
     Y="-3.355847048000E+06" />
   [..]
  </GCPList>
  <VRTRasterBand dataType="Byte" band="1">
  <Metadata />
   <ColorInterp>Gray</ColorInterp>
   <SimpleSource>
    <SourceFilename relativeToVRT="1">input.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="5409" RasterYSize="5408" DataType="Byte" BlockXSize="5409" BlockYSize="1" />
  <SrcRect xOff="0" yOff="0" xSize="5409" ySize="5408" />
  <DstRect xOff="0" yOff="0" xSize="5409" ySize="5408" />
  </SimpleSource>
  </VRTRasterBand>
  </VRTDataset>

Then I warp the input image to the final geotiff one by using

     gdalwarp -order 2 -co "TFW=YES" input.vrt output.tif

but the problem is that my output.tif has size (5419, 5414) instead of (5409, 5408), what am I doing wrong?


texts: Selecting features using an attribute which falls within a range

texts:
I want to select batches of features using an expression which gives a range of IDs to select.

Instead of writing this
"ID" IN (1,2,3,4,5) I would like to write something like this "ID" IN (1:5)

I tried the following but they didn't work:

"ID" IN (generate_series(1,5,1))

"ID" IN (array_to_string(generate_series(1,5,1)))


texts: creating many Shapefiles from from a single multi-polygonal shapefile using ArGIS tools

texts: I have a country shapefile with many districts polygons which are administrative boundaries, I want to have a shapefile for each of the district, which tool in ArcGIS can I use to accomplish this


texts: Spliting Curved Polygon with various length of lines using PyQGIS?

texts:
To split the curved polygon (river polygon) with various length of lines while keeping the original line position, I tried the procedure as bellow.

extend the lines using 'qgis:extendlines'

intersect the lines with polygon using 'qgis:intersection'

However, non-intended extended lines appeared compared to original lines (lines extending from a curved river), so I selected the extended lines based on original lines using 'native:extractbylocation' with 'overlap'.  However, no features were selected after the 'native:extractbylocation'

Can I get the features overlapped with original lines ?

This is my pyqgis script.

# hn_ln is original line layers
# hn_pg is polygon to be cut

output_line0 = processing.runAndLoadResults('qgis:extendlines', {'INPUT':hn_ln, 'START_DISTANCE':2000, 'END_DISTANCE':2000, 'OUTPUT':'memory:'})

output_line1 = processing.runAndLoadResults('qgis:intersection', {'INPUT':output_line0['OUTPUT'], 'OVERLAY':hn_pg,  'INPUT_FIELDS':[], 'OVERLAY_FIELDS':[],'OUTPUT':'memory:'})

output_line3 = processing.runAndLoadResults('native:extractbylocation', {'INPUT':output_line2['OUTPUT'], 'PREDICATE':['5'], 'INTERSECT':hn_ln, 'OUTPUT':'memory:'})


texts: QGiS 3.0: Line length without overlapping parts

texts:
I have overlapping lines in my layer and want a total length of the visible line without the overlapping part.

i.e.

I have two lines

First one goes A-B-C

Second one is B-C-D

How do I get the length from A-B-C-D?


texts: Extent and Projection on ASTER issue on my ArcGIS 10.1 ArcMap?

texts:
On my previous post on my ASTER issue, ASTER extent or projection issues. I decided to revisit on this issue. I am still puzzled as to why my ASTER has geographic coordinate system on its extent but the Spatial Reference tells me it is in UTM Zone 13 Northern Hemishphere and the datum shows me the Datum Not specified based on GRS 1980 spheriod. This is done in my office. The ASTER I download came from this website: USGS EarthExplorer

Now, I have ArcMap 10.1 at home and the extent and the projection shows me they are in the correct extent UTM numbers and correct projection.

How can one ArcMap 10.1 in my office and ArcMap 10.1 at home showed me two different extent numbers: One in geographic coordinate system and Projection UTM numbers ??

What could have cause them to show different extent numbers and why ? Could it be the installed change something within the ArcGIS program ?

For your information: My ERDAS 2011 show me the ASTER is in UTM meters and UTM Zone 13 projection and were able to overlay with polygons.

I will see if I can uninstall and reinstall my ArcGIS to see if it is the same results.


texts:
I have two exactly similar setups and I am trying to update my Asia layergroup with,
curl -v -u admin:geoserver -XPUT -d @Asia.xml -H "Content-type: text/xml" http://localhost:28080/geoserver/rest/workspaces/MapData/layergroups/Asia.xml.
In one setup the update gets done with no issues observed and I can see my modified Asia layergroup with new layers and their respective styles as defined in the Asia.xml.
But in the other setup I am getting "Invalid style group: No layer or layer group named 'Location' found in the catalog* Connection #0 to host localhost left intact" and the command does not update the layergroup as it does in the other setup.
There is no such layers or styles named Location in my setups and the same Asia.xml file is used in both setups.What could be the reason ?


texts:
I have a list of positions like this.

t, x, y, z

where t is a double representing time in seconds and x, y, z is a position expressed in ECEF coordinates.

I would like to interpolate the position between points in the list and I believe I can do that part, but given two points and time, what is the formula to do the interpolation?

lets assume the Position object contains time, x, y, and z

// @param p1 - first position
// @param p2 - second position
// @param time - time to get interpolated position for
private Position interpolateECEF (Position p1, Position p2, double time) {
   // need help here
}


texts: What's the theory (principle) behind shapely.ops.polygonize?

texts: I am learning Python package Shapely, but quite confused about the principle behind the polygonize function --- shapely.ops.polygonize. How does this function work to polygonize a set of lines (represented by the two endpoints)? What is the algorithm behind that?


texts: Attribute table reformatting fields on programme close

texts:
I'm creating a density map using the graduated style in qgis.  The data is created from the points in polygon tool, this is then saved as an integer (not sure why the plugin doesn't automatically do this!) and displayed as a style.  When the edits are saved and editing mode closed all the field types are retained, so far so good, however, when I close QGIS completely and come to open it again the fields all return back to string format preventing the graduated style.

Is this just a glitch?

How do I get it to retain the field type?


texts:
I have 10m spatial resolution image in my asset.
When I import into the window, perform kmeans segmentation and visualize the map, it looks like 100m resolution image.
How can I work with full resolution image?
Its important to see how the segments look at original scale.


texts: How to update selected rows of spatial data frame in R

texts:
I have a SpatialPointsDataFrame read via rgdal::readOGR

require(rgdal)
my_sp_layer <- readOGR(dsn="my_dsn", layer="my_layer_file")

# select rows using col_attrib column
my_sp_layer@data[my_sp_layer@data$col_attrib == 'some value',]

How can I update the values of the selected rows?

If it is a data.table, I would just use

my_data_table[col_attrib == 'some value', col_to_update := 'new value']


texts:
I'm a QGIS neophyte and I'm looking for a little help.  I'm trying to integrate QGIs into my forensic investigation practice, with regard to investigating and processing accident scenes.  I take a fairly high volume of photos and in some cases, it would be helpful to add a point layer to a QGIS map of the accident scene, which depicts where I took specific photos.

I'd like for the map to display the photo number, and align to the specific place the photo was taken from based on the photos EXIF data, including GPS coordinates, bearing, etc.,and to point, generally, in the direction the photo was taken, derived from the bearing.

The database I use to manage my images can separate out the EXIF positioning data.

What I'd love to be able to do is to have a Triangle or other marker that is positioned at the appropriate Lat and Long positions, facing the correct bearing, and labled with the appropriate photo number...

I doubt I'm the only person who would find this helpful.  Any ideas on where to look for explanations or tutorials on similar issues?


texts: Leaflet map shows no markers for GeoJSON

texts:
Following the example below I try to add an external GeoJSON to my Leaflet map, but its markers don't appear on the map. I use this file GeoJson. Any explanation?

https://maptimeboston.github.io/leaflet-intro/

  <title>A Cyprus map!</title>
  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.4/dist/leaflet.css" />

Error from console of Chrome:

$ is not defined

EDIT: The problem was the link of jquery and I didn't download any locally to.


texts: Converting GTiff to AAIGrid results in cell size switch?

texts:
I don't have a clue why the conversion of a tif to an asc has that kind of behavior:

original tif metadata

Size is 466, 427
Coordinate System is:
PROJCS["ETRS89_UTM_zone_32N",
...
Origin = (-1181164.876702783396468,4539238.589259737171233)
Pixel Size = (6486.620049172968720,7717.850831012333401)
...

It's hacky to convert it into a Ascii Grid because a VRT is needed

$ gdalwarp -srcnodata -1 -dstnodata -1 -of vrt grid_tif.tif junk.vrt
$ gdal_translate -of AAIGrid junk.vrt grid_asc.asc

The resulting Ascii Grid has that metadata

Size is 427, 466
Coordinate System is:
PROJCS["ETRS89_UTM_zone_32N",
...
Origin = (-1181164.876702783396468,7834760.894102003425360)
Pixel Size = (7075.195267921446430,-7075.195267921446430)
...

The original tif is created over gdal_grid of a point shapefile in UTM32N

gdal_grid -of GTiff -zfield value -l data_points -outsize 466 427 -a "average:nodata=-1" data_points.shp grid_tif.tif


texts: Mapping over a dictionary in Earth Engine?

texts:
Say I have a dictionary in Earth Engine as follows:

var dictionary = ee.Dictionary({
  key1: null,
  key2: null,
  key3: null
});

How do I use the map() function to change all the values of each key to a string "NoData" instead of null? According to the documentation for map(), I need to write a function consisting of two arguments: a key from a dictionary and a value, though I am unsure how to write this.


texts: Get a certain part of the project name in QGIS print composer label field?

texts: E.g. my project is named XXXX_YY_PRZZZ. Is there a way for me to get only the ZZZ part for my label field? Maybe with @project_basename or some variables like that?


texts:
I am using ArcGlobe.
I created contour lines using the contour tool, but some of them are not connected and are just fragments of lines, is there any tool that could connect these lines?


texts: Exporting a topo map with custom borders to STL

texts:
I'm looking to 3D print topographic models of all 67 counties of Pennsylvania, with each county individually printed. I am using the DEMto3D plugin according to this guide.

My issue is that I would like each print to be in the shape of the county rather than the rectangle that the source .tif file comes in. Any guidance you have to offer? My apologies if I can't get into the jargon of this process, I am quite the newbie.


texts: Getting a raster extent excluding nodata area

texts:
I have a GeoTiff raster in which half of its area is nodata. The problem is that the raster extent includes the nodata area.

Is there a way with GDal to recalculate the extent of the raster excluding the nodata areas?

UPDATE

Thanx to @user30184 here is how to get the raster extent of an image without nodata border.

Replace all nodata pixels with 0 and all other pixels with 1

gdal_calc.py -A source_image.tif --outfile out.tif --NoDataValue 0 --calc="1*(A>0)"
Get the shape file for the meaningful area.

gdal_polygonize.py -8 -b 1 -f "ESRI Shapefile" out.tif pols
Clip the source_image.tif with the shape file to create a new tif with updated extent.

gdalwarp -of GTiff -cutline pols/out.shp -crop_to_cutline -dstalpha 3559007_2014-05-08_RE5_3A_181809processed.tif OUTPUT.tif


texts: Faster way to get layer information in PyQGIS?

texts:
I know that you can quickly access the active layer with layer = iface.activeLayer()

Currently I'm working with multiple files in a Python script and my current method of accessing layer extents is this:

layer_file = "/path/to/file"
layer_meta = QFileInfo(layer_file)
layer_title = layer_meta.baseName()
layer = QgsRasterLayer(layer_file, layer_title)

layer.extent()

I'm just wondering if anyone knows a more efficient way to access the extents of a file without having to reproduce these four lines of code for each file


texts: Making world map to show where hole drilled through centre of Earth starts/ends on land vs water?

texts:
I am discussing drilling through the Earth with my students.

The classic cartoon version has a person drilling down in the United States and coming up in China, which has to be wrong because drilling from the northern hemisphere should come out in the southern hemisphere.

What I would like to get is a simple map of the Earth- single color background with continent / country borders only- with greyed-out / blued-out regions that don't have land on both sides of the dig.

I would like a map highlighting those regions of land that, when drilling through the Earth, end up on land on the other side.  I want to view/highlight regions of land/land, land/water, and water/water.


texts:
I am trying to embed a MapBox basemap that was created in the new studio (just came out a couple weeks ago) and when entering the MapID/URL and the Access Token in the CartoDB interface I get an error that says the URL is invalid. I am copying the URL that MapBox provides for sharing styles.

Do I need to pay for MapBox to be able to do this?


texts: Is it possible to export cross sections from QGIS to HEC-RAS

texts: Is it possible to create cross sections along a river in QGIS and then export the data to HEC-RAS?


texts: Choosing multiple options from QGIS Check Box or Value Map for feature attribute?

texts:
I am presently working with some layers containing location data for plots/parcel data. Some of the parcels have multiple use types.

Is it possible to configure the Check Box or Value Map (widgets type) so that it can allow me/user to input/choose multiple options in cases where the parcel has more than one use type?

I am using QGIS 2.18


texts: How to include colors in 16bits one band raster

texts:
I have many 16bits one band rasters (geotiff), who basically are in greyscale.
I had created a style for a color representation with "Singleband pseudocolor" in QGIS and save a qml file.
But for my customers, who don't all have geomatic services, it's easier to open a tiff with colors directly. But keep raster value is needed because some of customers want to use it.

Is it possible to include color in my tiff and if yes, how?

I can use QGIS, SAGA and python.


texts: Migrate datatypes from geometry to geography in MS SQL

texts:
I have three ESRI feature classes (POINT, LINE and POLYGON) that contain a geometry datatype column with data in a projected coordinate system. I'm looking for a way that I can convert & project from Geometry datatype to Geography.

Does anyone know of a tool or script I can use to convert all records permanently?

My research has show that I will not be able to do this with a simple SQL statement. So i'll probably need to use an Arc Toolbox or Python script or something.


texts: Prevent ArcGIS saving shapefiles into workspace

texts:
How can I stop ArcGIS creating output shapefiles as results of various tools?

I have got a script returning a certain value. However in order to get this value a have to use several functions such as arcpy.CreateThiessenPolygons_analysis or arcpy.Clip_analysis. These functions (tools) return maps and save them into my workspace. I would like to work with them in the script but I don't want them to be saved into the workspace since I won't need them any more.


texts: Can QGIS display individual WMS tiles as it fetches them?

texts:
I have a QGIS project accessing a base WMS layer with imagery.

How can I tell QGIS to display the individual tiles as it fetches them so that I won't wait for the whole image to load until I see something?


texts: Only show layer preview for group layers the user can access

texts:
I'm having problems using geoserver to serve geodata to some customers without displaying their respective group layers to each other (version 2.8.2).

Both users have access to their own workspace and it works well. The problem is that they can see a list of all group layers on the Layer Preview page. Even users that are not logged in at all can see the list. Can I prevent this somehow?


texts:
I have recently setup PostgreSQL and PostGIS on a debian server on a local network in my home, which is working perfectly.  I am able to access the db from other machines on the same network using the postgis connection in QGIS, which all works very well.

I would also like to access the same machine and postgres/postGIS database remotely from my office by making a connection through QGIS.  I have been struggling to understand the steps I need to go through to setup this kind of remote access.

Do I need a static IP for the home server?
Can I use dynamic DNS?
Do I need to setup openVPN? or similar?
I think I also need to forward ports on my router??

All the tutorials and guides I have found only really talk about the remote server being on the same network and not doing this over the internet. Any advice to pointing me towards some straight forward guides?


texts:
I'm pretty sure this is simple. I'm used to just doing a few maps at a time, and manually joining them doesn't take too much time. However, I've now managed to get myself into the situation where I have 20 infections I'm interested in. And 12 geographic areas. I'd like to create choropleth maps for each infection in each area. I don't think I can click manually through these.

I know that ArcGIS recommends this:

import arcpy

target_features = "C:/data/usa.gdb/states"
join_features = "C:/data/usa.gdb/cities"
out_feature_class = "C:/data/usa.gdb/states_cities"

arcpy.SpatialJoin_analysis(target_features, join_features, out_feature_class)

Is the standard "How to do spatial joins"
I've got that to work as a one off.

What I'd like to is iterate through the shapefiles in a file location (for each infection), join them to a shapefile (for the geographic area), save it as a new shapefile in a new location with a name like infection_area.shp

import arcpy,os

shpworkspace = r"S:\somestuff\maps\CSV Source file"
arcpy.env.workspace = shpworkspace
arcpy.env.overwriteOutput = True
shplist = arcpy.ListFiles("*.shp")
join_file = r"S:\somestuff\maps\areaname.shp"
pleth = r"S:\somestuff\maps\shpfile_areaname.shp"

#join shapefiles
for shpfile in shplist
arcpy.SpatialJoin_analysis(target_file, join_file, pleth)

I'm having a complete blank (as I'm new to the world of loops and python) on how to say "for each shpfile in shplist, do a spatial join to the join file, write it as pleth".

I've got a feeling I need to use os and split shpfile to create a name and location for the joined shapefile to be saved as.

I tried to build this in a model iterator but that nearly broke me.

I'm almost certain this has been replicated as a question, but how do I write a for loop that will do what I need?

(join a file in location a, to another in location b, write the output as partfile1name_partfile2name.shp to location b, repeat for all shapefiles in location a)

Solution:

import arcpy,os

#sets the workspace
shpworkspace = r"S:\somestuff\maps\CSV Source file"
arcpy.env.workspace = shpworkspace
arcpy.env.overwriteOutput = True
featureclasses = arcpy.ListFeatureClasses() #gets the feature class names from the workspace you set earlier

#the name of the template map
join_file = r"S:\somestuff\maps\areaname.shp"

# Join polygons to points
# This joins all .shp files in the workspace to the defined point file

for fc in featureclasses:
    name = fc.replace(".shp","_ward.shp") #create a name for the output
    arcpy.SpatialJoin_analysis(join_file,fc, name) #same join for each

I had to change the order to (join_file, fc,name)


texts: Colour matching satellite imagery from different time periods

texts:
Is it possible to colour match two separate satellite images taken at two different times?  I am working with two high resolution, multispectral images that were taken approximately one year apart and I want to do some change analysis with them.  Unfortunately one of the images is extremely washed out in urban areas and I'm trying to correct it (see attached samples).

Image From Time 1

Image From Time 2

Ideally, I would like to have both of the images to appear as similar as possible, preferably like the first one, so I can identify changes between the two time periods.  I am using ArcGIS 10.2 Professional, but if there are other options to deal with this I would like to hear about those as well.


texts: I couldn't find an option of "rule-based"-labeling in QGIS 2.10. For example that big cities are also labeled bigger than smaller cities. Am i missing something or doesn't it exist?


texts: How to add a geoserver layer to openlayers?

texts:
I am trying to add a geoserver layer with 2 points on the openlayers. But I only get the map without the points. This is what I try:

   <title>OpenLayers Example</title>


texts: Using in_memory in Python deletes values in field - Arcgis 10.4.1

texts:
I have a problem in using the in_memory function in a Python-script.

I have a shapefile test.shp with 3 fields (Alter, Geschlecht, Anzahl) and 2 rows of data. When I am using arcpy.CopyFeatures_management to copy the shapefile to in_memory it deletes the data of the rows.

I wrote a little test-script:

import arcpy

arcpy.env.workspace = "C:\\Python\\"

arcpy.CopyFeatures_management("test.shp","in_memory/test")
arcpy.CopyFeatures_management("in_memory/test", "test1.shp")
arcpy.CopyFeatures_management("test.shp","test2.shp")

You can download the shapefiles here: https://www.dropbox.com/sh/b47og9rjfevpi7v/AAB_YLz1l_eY3b-lFJPm9kP9a?dl=0

test1.shp is empty (rows are still here but without data) and test2.shp is an correct copy.

Can someone tell me what is going wrong?


texts: How to import KML file, with custom data, to Postgres/PostGIS database?

texts:
By custom, I don't mean a different structure of KML file, but in <description></description> section, I have data with HTML tags of this kind:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Document>
<Placemark>
  <name>New York City</name>
  <description><![CDATA[

Attributes

CITY_ID	150335

POPULATION	8244910

AREA	468.48

]]>
  </description>
  <Point>
    <coordinates>40.712964,-74.003886,0</coordinates>
  </Point>
</Placemark>
</Document>
</kml>

I want to easily import all this data, for each placemark, to a table with corresponding columns in the database, how to do this?


texts: Defining OID causes Query Layer on PostGIS to disappear from ArcMap?

texts:
I used query layer to load my layer from PostGIS and show it in ArcMap 10, I had two cases;

if I defined the OID in the configuration of connection, there was no problem with the attribute table of the layer but my layer disappeared from the map.
if I didn't define the OID, then the layer appeared in the map but when I try to open the attribute table of my layers the following error message occurs:

  Error reading OID form table. Reading rows has been stopped. Check
  that the datasource is valid. OID mapped column has null value

... in the table was empty.


texts: ArcGIS returns "IconSymbol3DLayer failed to load" when configuring a PictureMarkerSymbol

texts:
I'm trying to add a PictureMarkerSymbol to a SceneView, but I'm getting the following error upon reloading the page:

[esri.views.3d.layers.graphics.Graphics3DSymbolLayer] IconSymbol3DLayer failed to load (Request for icon resource failed:

Here is how I'm configuring the PictureMarkerSymbol:

var source = {
    geometry: {
        type: "point",
        longitude: 86.845,
        latitude: 27.991,
        z: 5194.164
    },
    symbol: {
        type: "picture-marker",
        url: "http://icons.iconarchive.com/icons/designbolts/free-multimedia/1024/Photo-icon.png",
        width: 20,
        height: 20
    }
};

var flayer = new FeatureLayer({
    source: [source],
    fields: [],
    objectIdField: "ObjectID",
    geometryType: "point"
});
map.add(flayer);

Can you please help me understand what am I doing wrong?


texts:
Given a simplified road network
with nodes (= intersections and turning points)
and uni-directional links (= path segments) between them,
I'd like to:

Know if there is a plugin or processing tool in QGIS

Identify / select the links crossing into a (closed) polygon (which outlines an enclosed area of interest)

Something like "terminal point is inside polygon"
NOTE:
I do not know the API of QGIS and want to save time, easily edit polygon and check the results or at least learn something, while avoiding to implement a solution using vector algebra.


texts: PostGIS/QGIS visualization with strange characters in attribute table

texts:
I loaded a shapefile into a PostGIS database. When I visualized it in QGIS I noticed that in the attribute table, some fields with integer data type, are displayed in a strange form.

However when I save the data as shapefile the fields look normal.

Why does this happen?


texts:
I have downloaded bathymetric data into QGIS. In order to use it I need to conduct a float to raster conversion. I have not been able to find how to do this on any of the forums. I have unsuccessfully tried using 'Raster' --> 'Conversion' --> 'Translate' option but still have found no attributes in my attributes table.

Here is a shot of the float file (grey shades file in the background) and the habitat area in green. I have also opened the attribute table of the habitat area file for you to see, which shows the all the statistical areas included, the sum of which make up the entire habitat area. Eventually i am planning to use interpolation plugin (TIN) to calculate HABITAT volume. That is the purpose of this exercise.


texts: Making raster layers calculated from point feature congruent?

texts:
Raster A and Raster B are both calculated from the same point feature. Raster A results from Natural Neighbor Interpolation in the Spatial Analyst toolbox. Raster B results from Point Density in the Spatial Analyst toolbox. However, the rasters are not congruent with each other:

Raster B has 2 more columns and two more rows than Raster A
Cells of Raster B are shifted relative to cells of Raster A

How to start in order to get congruent rasters from the same point feature? Shift of one of the rasters is not a solution as the value within a cell should be accurate, i.e. it should represent the area covered by the cell.


texts: Capture and parse output of lascheck from lidR

texts: To aid in quality control of a set of LAS files, I would like to write out all errors/warnings from lascheck(), for each LAS file (i.e. so not the LAScatalog implementation of lascheck(()). I'm looking for any tips to extract lascheck() warnings while looping through the entire list.


texts: Missing Layout and Data Frame Buttons in ArcGIS Desktop

texts: I have lost the Layout and Data Frame Buttons at the bottom of my map window in ArcGIS Desktop 10.3.1. How can I get them back?


texts: Plotting Dataset of netCDF using matplotlib.pyplot displays rotated plot

texts:
I have read a Dataset in NC file using netCDF4 in python and then tried to plot the dataset using matplotlib.pyplot. But issue is the plot shows the image but with 90 Degrees clockwise rotation. Following code I have used:

from netCDF4 import Dateset
import matplotlib.pyplot as plt

f = Dataset(r'F:\data\data.nc', 'r')
data = f['et'][:]

plt.imshow(data)
plt.show()

Kindly let me know what can be the issue.

Also one doubt I am having, I have 2 NC files having same dataset and in that, dataset in one file is getting properly displayed by using np.rot90(data) which was suggested to me in the answer given below and in another file I need not to do "np.rot90(data)" but instead I need to do "data[::-1]" to display the data correctly.... And interesting thing is when I am importing these 2 files directly in the Arcmap both are getting properly displayed....can u suggest what can be the issue? Whether it's the issue with NC data?


texts: Calculating Albedo for Landsat 8 using QGIS with GRASS?

texts:
I am trying to estimate albedo for Landsat 8 OLI. The idea was to use the GRASS function i.albedo within the processing toolbox of QGIS. However, the option for Landsat 8 OLI is not available, therefore I tried to use directly GRASS GIS. Being unfamiliar with GRASS I could not figure this out.

Steps I followed.

More Downloaded DN value images from UGSG Earthexplore website. (I want to use the Surface Reflectance Higher-Level Data Products from USGS which makes steps 2-3 redundant, but the values were off when I tried. I had read that the values stored in this product needs to be multipled by a factor, ex. 0.0001, so that I can use this product?)
Performed TOA (Top-ofAtmosphere) reflectance conversion from DNs (wanted to use GRASS function i.landsat.toar but could not figure this out, used semi-automatic classification plugin for QGIS instead.
Wanted surface reflectance so I applied DOS1 correction within the above plugin.
Cannot calculate from surface reflectance to Albedo for Landsat 8?

If I could do this within QGIS it would be best, or if someone is familiar with GRASS GIS, I could use a simple tutorial.

Remote Sensing Landsat Surface Reflectance and Albedo


texts:
I'm building an interactive map using CartoDB.js in an HTML page. I have the user input an address that gets geocoded. I then use the user's lat/lon coordinates to perform SQL queries on CartoDB. The first query looks for a polygon that intersects that lat/lon, then selects three other points that share and ID number with that polygon. It looks like this:

WITH tstations AS
    (SELECT * FROM dsny_refuse_recycl_net_wgs84
    WHERE districtco =
        (SELECT borocd FROM nycd
        WHERE ST_Contains(the_geom,  ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326))
        )
    )
SELECT * FROM tstations

This returns me the three points that I want, as "tstations".

What I wanted to do next is create three lines, connecting the user's lat/lon point to each of the three tstations points, and to display these lines on the map.

I'm using the following query, which works on the CartoDB website's SQL window:

WITH tstations AS
        (SELECT * FROM dsny_refuse_recycl_net_wgs84
        WHERE districtco =
            (SELECT borocd FROM nycd
            WHERE ST_Contains(the_geom, ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326))
            )
        )
SELECT null as cartodb_id,
ST_MakeLine
   (ST_Transform
       (ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326),3857),the_geom_webmercator)
       as the_geom_webmercator
FROM tstations

(I'm basing myself off this tutorial: http://docs.cartodb.com/tutorials/query_by_distance.html)

When I use the first query above within my HTML with CartoDB.js, it does work, and returns me the three points I want. But when I use the ST_MakeLine query, it doesn't work - CartoDB's API returns a 400-Bad Request error. This is what it looks like inside my HTML script:

function updateQuery() {
    sublayers[1].set({
      sql: "WITH tstations AS (SELECT * FROM dsny_refuse_recycl_net_wgs84 WHERE districtco = (SELECT borocd FROM nycd WHERE ST_Contains(the_geom, ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326)) ) ) SELECT null as cartodb_id, ST_MakeLine (ST_Transform (ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326),3857),the_geom_webmercator) as the_geom_webmercator FROM tstations",
      cartocss: "#dsny_refuse_recycl_net_wgs84_copy_copy{ line-color: #FF6600; line-width: 2; line-opacity: 0.7; }"
    });
}

Does anyone know what I'm doing wrong here? I assume it has something to do with using sublayers.set to generate the new lines, but I'm not sure.

Here's the API error response:

{"errors":["Postgis Plugin: ERROR:  syntax error at or near \"AS\"\nLINE 1: ...(ST_MakePoint(-74.0059413,40.7127837),4326),3857) AS the_geo...\n                                                             ^\n\nin executeQuery Full sql was: 'SELECT * FROM (SELECT d.cartodb_id, ST_MakeLine( d.the_geom_webmercator, ST_Transform(ST_SetSRID(ST_MakePoint(-74.0059413,40.7127837),4326),3857) AS the_geom_webmercator FROM dsny_refuse_recycl_net_wgs84 AS d JOIN borocd AS b ON d.districtco = b.borocd WHERE ST_Contains(b.the_geom, ST_SetSRID(ST_MakePoint(-74.0059413,40.7127837),4326))) as cdbq LIMIT 0'\n  encountered during parsing of layer 'layer1' in Layer"]}


texts:
I'm working on an application where a desired functionality is for users to start with the big picture (i.e. a single map layer), and gradually add more and more conditions to display only the area that satisfies all of them (i.e. the intersection of all the selected layers). I've been able to do this with Javascript on the client, but this has some performance issues and ideally I'd like to process this on the server.

The data I'm working with has already been prepared as polygonal layer or raster layer mapservices. Can I pull data directly from these services for processing?


texts: Having issues getting arcpy.Project_management to execute

texts:
Python keeps returning me with a Error 000208 on the Project (Project_management()) function. I've checked for improper syntax for all arguments but nothing is striking me as wrong.

###Import modules into the script
import sys,arcpy,traceback,os

from arcpy import env
###Set environments
inwrkspc = env.workspace = r"E:\WWU_2017\ENVS423\Assignment_9\Original_Data\Ex9_Data\Ex9_Data_052516\InputProjections.gdb"
outwrkspc = r"E:\WWU_2017\ENVS423\Assignment_9\Original_Data\Ex9_Data\Ex9_Data_052516\OutputProjections.gdb"
print "Workspaces Set"

###For a feature class item in InputProjection.gdb...
fcList = arcpy.ListFeatureClasses()
for infc in fcList:
    print infc
    descSRN = arcpy.Describe(infc).spatialReference.name
    print descSRN
    ###if the descSRN(Name) and output proj string is true then pass/ignore this file
    if descSRN == "NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet":
        print "This data are already projected in the correct coordinate system."
        pass

    ###Otherwise, run reproject
    else:
        print infc
        outfc = os.path.join(outwrkspc, infc)
        ###Output coordinate system
        outcoordsys = arcpy.SpatialReference('NAD 1983 HARN StatePlane Washington South FIPS 4602 (US Feet)')
        ###Project tool exported from ArcGIS
        arcpy.Project_management(infc, outfc+"_06052017reproj", outcoordsys)
        print "{0} was reprojected to NAD 1983 HARN StatePlane Washington South FIPS 4602 (US Feet)".format(infc)

ArcPy ERRORS:
ERROR 000208: Error creating output feature class
Failed to execute (Project).


texts:
I would like to use QGIS to create my map. One data feed is the gtfs regarding to bus routes. I can view it on Google Earth (although laggy), but when I tried to import the .kml into QGIS, it just ran out of RAM (for 30 min) and crashed my map. If I import 10 in the 18xx layers, the route is duplicated 10 times (Google Earth could show all the routes, although laggy).

My procedure

I am trying a new method, so...

What I did was download the .txt on TfNSW website (Google tfnsw data and login), then use a .py script called Transitfeed by bdferris to convert the .txt into .kml. Finally import the .kml as vector layers, but I did viewed in Google Earth prior that so I though it will be okay.

I am intended to convert the shape.txt into route shapes.

Problem

I was thinking of the possibilities and solutions, for example:

Faulty .kml
Faulty .txt in the gtfs
Faulty conversion from .txt to .kml

I was thinking if either first two possibilities are valid, then what is the suitable validation (of my .kml) solution? Otherwise I would like a cost-less (i.e. monetary) solution for viewing gtfs .txt file into QGIS.


texts: Using Geoserver extension CSVDataStore alters existing data store types

texts:
I am trying to create a custom extension into Geoserver which will act as a datastore.

Being new to GeoServer, I have started with the CSVDataStore tutorial from the GeoTools website.  I have managed to compile the extension and installed it into GeoServer.

The problem is that as soon as I create a new data source based on the new extension, the existing datastores based upon shapefiles change their type to CSV as well.

Anyone experienced this before?

Edit:  I have switched from version 2.11.1 of GeoServer to 2.10.3, and so far I have not seen this behavior. I am not sure yet if this actually is the solution, but I will keep this post updated


texts: Leaflet gives javascript error after rendering

texts:
I get javascript errors but the map is visible without any problems. My code
ts file
options = {
layers: [
  tileLayer("https://tiles.stadiamaps.com/tiles/alidade_smooth/{z}/{x}/{y}{r}.png",{
    maxZoom: 2,
    minZoom: 2,
    noWrap: true,
  })
],
zoom: 2,
zoomControl: false,
center: latLng(49.3824, 18.1055),
};

HTML

Any ideas why?


texts:
I want to display an image on the feature popup in qgis2web. I use the latest version of qgis2web and Qgis 3.0.3.

I chose Attachment as the widget type and added the path to the photo using the attribute table. But in the qgis2web export only the path to the photo displays and not the photo itself. It worked in Qgis 2.x using the widget type Photo.


texts:
I just edited a shapefile in ArcMap and I was unaware that this would edit the shapefile in ALL ArcMap documents in which I had used it.

I want to restore this shapefile but cannot figure out how.

I tried on Windows "view previous versions" but it says there are no previous versions.

I also tried according to the ArcGIS website: http://desktop.arcgis.com/en/arcmap/latest/extensions/production-mapping/restoring-a-previous-version-of-a-file.htm
But I cannot manage to get the product library window the way they suggest.

What can I do?


texts: Problems with cleaning the topology of a shapefile using GRASS tools and snapping tools

texts:
I have a shapefile with almost 1000 polygons (each polygon representing a farm). There are a lot of topological errors (gaps/slivers and overlaps between polygons). Below is a "zoomed in" example. I am using QGIS 3.10.

I tried out the GRASS "v.clean" function as well as the snapping tool in QGIS and have two questions:

Using "v.clean", which specific tools should I use to remove the gaps / overlaps and how can I determine the threshold value that I should use (where applicable)? I tried using "break", "rmdupl" and "rmarea" (using various different threshold values), but I do not see any changes in the cleaned output (the gaps / overlaps remain). What reasons could there be that this is not working?

Using the snapping tool, how can I move an entire selected polygon (yellow on the map) to match up with the borders of the adjoining polygons, instead of snapping each and every vertex to each other? It looks like the entire selected polygon just needs to be moved, and the gaps and overlaps will be eliminated automatically. Furthermore, I am worried that if I start snapping the segments/vertices, then the shape of the polygon will change if I do it incorrectly (and it will take a lot of time doing this for the entire farm shapefile!).

I am investigating the movement of animals in relation to farm boundaries, so it is vital that the boundaries of the polygons do not change shape.
Here is a link to my farm shapefile
Please note that I am a QGIS beginner, so I need an explanation that is not "too advanced".


texts: Imported CSV data is not displayed in QGIS

texts:
I watched a simple video tutorial on how to accomplish it, but nothing shows up on my project when I add the layer. Even if I try zoom to layer, the screen just shows blue. I'm not sure if it has something to do with my Excel, computer, or program. I really need help as I am new to QGIS.
(I just posted this, but it won't let me comment to people under that question any more)
I have a layer of google hybrid, but when I do add the new layer it is at the top and checked transparency and everything. Its just not accepting the data I entered into the excel spreadsheet is my guess to what the issue is.

I have the spreadsheet with the name, longitude, latitude.
Save as CSV
In QGIS go to Layer > Add layer > Add Delimited Text Layer
Upload the CSV Excel file
Add

Nothing shows even tho there is a line in the layers (bottom left) showing that "something" is there but if I choose "zoom to layer" it changes to a blue screen I have a
Base layer Google Maps hybrid, and an area of interest on Georgia


texts: How to pull data from multiple databases via Qt .ui form?

texts:
Using a single .ui form and wthin QGIS, can I show multiple records from another table that is associated with a single polygon?

Table 1 contains planning applications ref
Table 2 contains multiple address associated with one planning ref.

I click on a polygon that has a planning application reference and within the same form, I want to view all addresses associated with it?

I am pulling the information from a non-spatial oracle database and my polygons exist in esri shapefiles.


texts:
I am trying to display GRD and GRD.gi files in QGIS with no luck. The files seem to load fine using Mapinfo, but mapinfo doesnt allow you to convert to geotiff. Any ideas?
This is a geophysics data set from a government website dating back to the 1990's


texts:
I have a shapefile in QGIS 3.16 that I created by performing a spatial join.  It displayed properly on the map. I then edited some of the data in one of the spatially joined fields. After doing so, some of the features don't show up on the map. It seems random, though. Some unedited ones are visible, some are not. And vice versa.
I've searched online for an answer and have not found anything. Might it have to do with my spatial joins layers not having a spatial index? Or I also noticed the spatial join layer did not have a CRS.
I'm lost.


texts: Find out which Python executable QGIS is using on Mac OS

texts:
I'm trying to find the path to the Python executable QGIS 3.6 is using on Mac OS. Usually (i.e. Python "outside" of QGIS) I'd:

import sys
print(sys.executable)

But that only gives me the path to QGIS:

/Applications/QGIS3.6.app/Contents/MacOS/QGIS

Is there any straightforward way to find this out?


texts: Can anyone recommend a GUI for arcpy? So if I want to run my script as an exe with some blank text boxes(or text boxes with default entries), radio buttons, check boxes etc.


texts: Converting geometry to WKT format by using QGIS Expression

texts:
I have a gpkg file and I am changing the geometry to wkt format by using expression geom_to_wkt( $geometry ) in the Field Calculator.
The output looks like this in the image. The geometry type is not written in capital letters and the polygon is not closed as WKT format is defined here http://giswiki.org/wiki/Datei:Geometrietypen_im_WKT-Format.jpg.
Am I doing something wrong?


texts: Line and polygon intersection in GeoPandas/Python

texts:
I am trying to intersect two GeoPandas data frames, the first one is the street data (streets) and the second one is a buffer (polygon) around a point (buffer).
I want to find the intersection of streets in this buffer zone and if a street extends outside of the buffer I want it to be cut on the border. I can do the following:

streets_in_buffer=geopandas.sjoin(streets,buffer,how='inner', op='intersects')

However, I cannot cut the street (in red) at the border of the buffer (in blue) and it extends outside of the boundary like this:

I was wondering if there is a way to only get the portion of the streets inside the circle and end up with something like this:


texts: Editing vector data over the net: WFS or PostGIS?

texts:
I'm using QGis to update the attributes of a large vector dataset (~20,000 polygons) stored in PostGIS. I'm going to have to do that to other large datasets so the idea is to bring in more people, each editing a separate table (no concurrent editing).

I'm wondering what would be the best option: Have the users connect directly to PostGIS and edit the data? Or publish the data in GeoServer and use WFS? Are there any gains in speed if I use WFS (cache, load only polygons in visible window)? Or would the extra layer (geoserver) slow me down? I'm using QGis as a desktop client.


texts:
As part of a project for a wind turbine assessment, I would like to make a wire frame image to compliment my Viewshed / line of sight and Photo montage.

I have never created a wire frame image from dem data and haven´t got a clue where to start.

Could anyone direct me in the right direction regarding software, preferably OpenSource?


texts: How to identify each circles in Google Maps API?

texts:
By the Google Maps Drawing API, I can create lot of circle and other shapes. addListner method helps to get the radius, latitude and longitude of the circle. Now I want to identify each circle for storing its coordinates and radius. The intention is when user clicks in circles, all these values stored in the input button for each circles and submit to database.

How can I move forward ?


texts: TopoGeometry commands with postgres are really slow

texts:
After messing around quite a while, I was finally able to get a topogeometry converted from my geometry table following the example of toTopoGeom. I have the topology extension, with a topology schema (SRID = 5367) and a table in the public schema with a topogeometry column.

However, after running a couple of minutes I quit the query because it would take too long. Now I have around 4000 rows in the topogeometry table and only ~10% converted.

Then I wanted to convert the topogeometry to a topojson following the example of the AsTopoJSON. It works ok, but it takes 7 seconds to generate 4000 rows, which is really slow in my opinion.

So the question is; is this a common problem and is there an alternative approach to speed things up?

used code

Insert topogeometry from geometry:

INSERT INTO grid_topo(topo)
SELECT topology.toTopoGeom(geom, 'grid_topo', 1)
FROM grid
WHERE gid >= 0 AND id < 500;

(added the condition to limit the running time).

Generate TopoJSON from topogeometry:

SELECT AsTopoJSON(topo, null) from grid_topo


texts:
I update from ArcMap 10.3 to 10.4.

When I export a map, I always used the "normal" north arrow, i.e. the first arrow on the top-left in the image shown here:

After I update to ArcMap 10.4 is that now, when I select the first north arrow type, that in the preview still appears as the one I always used, something different pops up on the map. The new north arrow looks a lot like the central one on the figure I showed before.

Am I the only one with this problem? How can I go back to my previous north arrow (with just the "N" letter written on it)?


texts: Advanced Digitizing Panel bad value or bug?

texts:
This is my first post here.

What would cause the "Advanced digitizing panel" D (length) to show up wrong?

I'm in the US using WGS 84 as CRS.

If I add a feature (line) Advanced Digitizing Panels reads "d" value says 10.6, but using the measurement tool is 24ft and confirmed correct with Google Earth its 24ft.

Ultimately, my goal to know the length the line I'm adding.

Any suggestions?


texts: Showing shorter list from long domain list?

texts: I have a geodatabase with a feature class that has several domains attached to various attributes in the table. A couple of the attribute have very long domains associated with them, where much of the data is not used any more. For example: one of the fields is for the data observer; we currently only need four of the people listed and not the entire list of all the data collectors for the past 12 years (which is like 30 people). It would be much more convenient to just be able to have a shorter drop down list for the field data collectors.


texts:
I'm trying to install QGIS 3.x on Raspbian Stretch with Desktop - I successfully installed QGIS 2.x using the QGIS sources as described in their install instructions.

I found this build-from-source guide using this forum post, but my linux-fu is weak. I'm wondering if anyone has found a simpler way of getting a standard QGIS 3.x desktop install on Raspbian Stretch?

I'm particularly confused by the 'build Debian package' option when building from source - what previous sections do I skip (if any)?


texts: Download openstreetmap data filtered by specific changeset comment

texts:
From openstreetmap, I want to download only the features that have specific change set comments (like #ABC). I already know how to download OSM data (using JOSM) and found this webpage shows filtered result only; no download option available (my required format) .
It would be a great help for a newbie if anyone provides solution.


texts: Creating edge lines between specific polygons

texts:
I have a polygon layer containing habitat types. It looks something like this:

where dark green is forest, yellow are open lands and the other signatures different other habitat types. In order to get something like a forest edge-line, I'd like to create a line shapefile representing those linear structures where forest and open land touch each other. I could of course use the polygons to lines tool in QGIS and then select the forest edge lines by hand - however, my data are too large for this procedure.

Any suggestions on how to approach this problem in QGIS or R? Please note that - unfortunately - I am not familiar with python at all.


texts: Choosing UTM zone to use for large country?

texts:
I want to reproject OpenStreetMap roads data from the current projection (WGS 84) to UTM, since I read in another GIS SE Q&A(Getting $length in meters in QGIS?) that WGS 84 cannot be used to measure length in m or km.  My end goal is to calculate the total roadway kilometres within each district, so that I can enter that variable into my regression model.  The country in question is Indonesia, and as you can see from the map below, it covers UTM zones 46-54.

Image source: Wikimedia, http://upload.wikimedia.org/wikipedia/commons/e/ed/Utm-zones.jpg

Does it really matter which one I choose?  Can I choose the one that's roughly in the middle and call it good (that would be zone 50) or should I pick the zone that is more densely populated with roads (Jakarta/Java, zone 48)?

I am only barely functionally literate in GIS.  The only formal training I had was a 3-day course back in 2008.  Really, the only reason I need the roadway kms variable was because one of my thesis committee members wanted me to account for it in my model...

And to confirm, once I have reprojected to the UTM CRS, I can use my GIS software's built-in length tool to directly calculate the length in metres (or divide by 1000 for km), right?

I am using QGIS.


texts:
I have a line layer that is rulebased labeled to show green segments and red segments based on 'status' attribute.
Currently, I download data from an excel sheet on sharepoint as a csv and then import it to qgis and create a join. Finally I'll save the new shapefile in postgis and publish on qgiscloud.
Is there a way to scrape the excel data on sharepoint and have qgis cloud prject the changes live?


texts: Display of DIVA-GIS Gazetteer *.dbf Files in QGIS

texts:
I've looked through other questions/answers and not seen what I am looking for...

How do I display the contents of a DIVA-GIS Gazetteer *.dbf file in QGIS?

Ideally I would like to take contents of the dbf and display a symbol according to coordinate(s), defining symbol character from another column and finally labelling according to a fourth column.


texts: Getting top 1 record from shapefile by descending order using ArcPy SearchCursor

texts:
I have a shapefile from the date field in it. I'm looking to get top 1 record date value by descending order.

readRows =    arcpy.SearchCursor(shapefile, "", "", "Last_Updated" ,"Descending")

Is descending order possible from ArcPy?


texts: Draw lines and Move the marker at the same time with KML format

texts:
it is possible to draw a line and put the Marker in the top of the line with KML format? or have i to make two KML files, it's that to say, One KML file for the Lines and another one for The Markers?
We are about to develop a GIS track app, and we are testing what is the best way to achieve our project.

here is and example, about what we want to do:
GIS Track Example, Moving line and Marker

Well is good to explain, that we gonna get a new KML or the data every 20 seconds to redraw the lines and the markers in the map.

Or do you have better ideas? KML, Database(Mysql), Xml?

Thanks to all.


texts: Projecting MODIS LST (MOD11_L2) Data using GDAL (+ Python)

texts:
I am trying to analyse some MODIS LST data (MOD11_L2) and am having some difficulty projecting the files into WGS84. I know that I can use the MODIS Reprojection Swath Tools, but I am trying to do the whole process using gdal (which I am currently running through python as this is where the rest of my process is).

In python, I have been using:

import subprocess
subprocess.call(some gdal functions)

I have tried the gdal_warp function but am having trouble with the parameters of the tool. Here is an example file:
http://e4ftl01.cr.usgs.gov//MODIS_Dailies_F/MOLT/MOD11_L2.041/2015.08.15/MOD11_L2.A2015227.0215.041.2015228093331.hdf.

I have tried looking at the help here: How to reproject MODIS Swath data to WGS84 but it does not seem to work with the LST data which I am using.

I want this to work so I end up with a Geotiff file in WGS84.

I'm fairly new to both gdal and python!


texts: ArcGIS Geoprocessing Service Output Location

texts:
I am publishing an ArcGIS geoprocessing service using a Python script tool and can't seem to figure out the correct output path to use. I would like to share the results as a map service, which means that, according to this documentation, I should not use the "in_memory/output_file" scheme.

So, my script writes to a permanent location on my server:

output_path = r"\\server\path\to\data\geodatabase.gdb"
output_layer = os.path.join(output_path, "Result")
arcpy.SetParameterAsText(1, output_layer)

But there are some issues that seem to remain.

Although the tool runs successfully in ArcMap, it fails after its published on the server with an error that says the output layer already exists.
Assuming I set the geoprocessing environments to overwrite the output, is this the behavior I want? What if there are multiple instances of the geoprocessing tool running... will they not conflict with each other when trying to read/write?

If I understand correctly, a new folder is created for each job, something like:

c:\arcgisserver\directories\arcgisjobs\%service%_gpserver\%jobid%\

I just can't figure out how I can define my output file to be placed there so that it can get pushed up in a temporary map service result.


texts: GDALWarp and Target Aligned Pixels Not Effective

texts:
I have two sets of TIFFs of different resolutions (100m and 1000m),  different alignment, and varying data types with corresponding no data values.  I need to downsample the 100m TIFFs to match the resolution and alignment of the 1000m TIFFs.  I also need to preserve no data values in each TIFF.

I am having a heck of a time accomplishing all three objectives in a straight forward manner without using multiple tools with repetitive resampling degrading the information.  Ideally, I would like to be able to accomplish all three objectives from the command line via batch, stand alone Python, or the Python window.

I have been able to achieve the correct cell size, no data value, and alignment, but by using awkward combinations of

GDAL_Translate

From inside QGIS and from a batch file running
c:\OSGeo4W64\bin\GDAL_Translate.exe;

Align Rasters

Only from inside QGIS where the menu does not allow the specification of no data values for both the source and reference TIFFs leading to edge-effects as 100m TIFF no data values are included in averaging during resampling.  Will work only if 100m raster has already been resampled to 1000m and is then resampled again by Align Rasters.  I have not been able to run Align Rasters from any sort of command line in batch or Python.  I have not even found Align Rasters on my computer as Python nor as executable.

GDALWarp

Documentation seems to suggest that I could use GDALWarp specifying Target Aligned Pixels (-tap) to accomplish all three in one go.  I can run
c:\OSGeo4W64\bin\gdalwarp.exe
from a simple batch file, but the -tap parameter is ineffective at aligning the result.

What I have tried is,

C:\OSGeo4W64\bin\GDALWarp -tr 1000 -1000 -tap ..\ROMS\BottomSalinity_max.tif -r average -srcnodata -32768 -dstnodata -32768 -of GTiff -overwrite bathy.tif bath_1000m_ROMS.tif

Can anyone see what I would need to change here?

I have to confess I am not sure I properly understand the single sentence in the documentation describing the purpose of the -tap parameter.

"-tap:  (GDAL >= 1.8.0) (target aligned pixels) align the coordinates of the extent of the output file to the values of the -tr, such that the aligned extent includes the minimum extent."

Other threads (and I have read a lot of them recently) suggest that -tap is the way to align rasters.  Am I just failing to supply the correct parameters or, is -tap not the way to align rastes?

Is there a way to access Align Rasters from the command line avoiding gdalwarp?  I do not need to warp -- just downsample and align.  I am only looking at gdalwarp as a tool through which I can successfully set the parameters for downsampling and aligning.


texts: Including libfcgi files when installing QGIS on SUSE Linux

texts:
I'm using SUSE Linux Enterprise Desktop 11 SP2. When I try to install QGIS from YAST, it seems that the package "libfcgi.so.0" and "libfcgi.so.1" are missing.

Where can I get those files?


texts: Relationship between ST_SubDivide and ST_Union?

texts:
From the docs on ST_SubDivide()

  Turns a single geometry into a set in which each element has fewer than the maximum allowed number of vertices. Useful for converting excessively large polygons and other objects into small portions that fit within the database page size. Uses the same envelope clipping as ST_ClipByBox2D does, recursively subdividing the input geometry until all portions have less than the maximum vertex count. Minimum vertice count allowed is 8 and if you try to specify lower than 8, it will throw an error.

If one does this..

SELECT ST_Union(geom)
FROM (
  SELECT ST_SubDivide(geom)
  FROM tbl
  WHERE id = 1
) AS t;

Should the resulting ST_Union( ST_SubDivide(geom) ) be the same? What exactly does envelope clipping mean? Are we losing area or detail?


texts: Is it possible, in QGIS, to create buffers with a set radius, but with individual values taken from the attribute table?

texts:
I am looking to create an output that consists of an 11KM buffer around a number of discrete  points, where the buffer has a value defined by a column in the data sets attribute table.  Is this possible?  If so, is it possible to aggregate the buffer values that 'overlap'? The points represents the origin of each buffer, where the buffer has an 11Km radius.  I would like to assign each circle (I have made the fill 95% transparent to simulate, visually, what I am looking for) a value equal to a column in the attribute table, and then aggregate the values wherever two or more buffers overlap/  So, for instance, the two buffers on the souther edge of the city would have an aggregate value equal to the sum of each of those buffers.  The buffers in the middle of the city would have an aggregate value equal to the sum of all overlapping values. In essence, the darker the shade, the higher the value.

Does this make sense?

BTW, thanks @underdark for the suggestion of adding an image.

I was able to add a buffer to the points (as described below) and then by joining the buffer values to the attribute table I was able to define a value for each buffer and visualize, as seen in the following image.  Now, I am uncertain how to aggregate the values of each buffer with those that overlap it. I suspect that I need to create a raster of values, but am uncertain.  Could anyone please suggest a method?

Thank you!


texts:
I'm trying to add a geoJSON layer to my OpenLayers3 map following this example.  Below is my GeoJSON

{"type":"FeatureCollection",
 "crs":{"type":"name","properties":{"name":"EPSG:4326"}},
 "features":[
  {"type":"Feature",
    "geometry":{
      "type":"Point",
      "coordinates":[-78.6569,37.4316]},
    "properties":{
      "Title":"Test ",
      "Assessment":"Test",
      "State":"Virginia",
      "Brief Description":"This ia a brief Descript",
      "Details":"

Here are the details

"}}]}

When I add this layer to the map I get a point at [0,0] instead of the coordinates in the json file.

function loadGeoJSON(geo)
{
    console.log('loading Geo JSON');
    console.log(JSON.stringify(geo));
    var vectorSource = new ol.source.Vector({
        features:(new ol.format.GeoJSON()).readFeatures(geo)
    });
    var vectorLayer = new ol.layer.Vector({
        source:vectorSource,
        style:styleFunction
    })
    map.addLayer(vectorLayer);
    for(var i = 0; i < vectorLayer.getSource().getFeatures().length; i++){
        console.log(vectorLayer.getSource().getFeatures()[i].getGeometry());
    }
}


texts: Image collection images to a 1D vector where there is one value corresponding to each image

texts:
I am currently working with MODI's land surface temperature image collection dataset of 20 years of data. I want to convert each image to a single value representation using mean/max function. I tried the map function, but it returns only images, whereas I am looking for a vector of float values where each value is the max/mean of the temperature of that region.

How can I solve this issue?


texts:
Using ogr2ogr command line I can overwrite or append to an existing PostGIS table like:

ogr2ogr -f "PostgreSQL"-overwrite PG:"host=myserver port=54321 user=user dbname=mydbname password=mypassword" C:/Temp/my.shp

Next, I would like to incorporate this logic into a python script.  I've found some useful sites, however I'm not finding any examples using these methods.

Open Source RS/GIS Python Week 1

PostGIS OGR Examples

Append a shapefile to a postgis table using the GDAL/OGR Python interface

How do you overwrite or append to a PostGIS table using ogr2ogr commands in python?


texts: Clipping/Splitting line layer based on another layer line running in parallel/alongside it using QGIS

texts:
I know how to Clip layers files based on another layers thanks to the Geoprocessing toolbox but how do I clip/split a line layer based on another line that is following parrelel/allong side it?
Look at the picture below, the Blue line is seperate layer following along the green on in certain points, I want to "split" this blue line layer where it follows with the green one.
So that only the sections where the green layer line does not go remain (highlighter in Yellow).

How do I do this without manually editing?
It is a process that I want to automate into a model.
Currently working in QGIS 3.4.5.


texts:
I have a shapefile with 10,000 polygon features in QGIS I want each polygon to be split into multiple polygons features (partitioning whole polygon area) using some attribute value say NoOfPartions. I am using QGIS for same but there doesn't seem any way to do so for all features. Manually clipping is not feasible due to large number of features.

Is there any possible way to do so in QGIS/PyQgis?


texts: How to convert SRTM's Lat/Lon to Landsat's utm reference system in GDAL?

texts: SRTM data is given in Lat/Lon coordinates, while Landsat is given in UTM coordinates. Is there a way to reproject the SRTM data into the Landsat CRS in GDAL (Python)?


texts: Creating Buffer of specific area using QGIS?

texts: I am trying to create a buffer which would include only the Mediterranean coast of Spain (excluding borders), running from the south western extremity close to Portugal, to the north eastern extremity close to France. Is there a way to to this based on the map I have?


texts: Difference in raster pixels between Windows x86 and Windows x64

texts:
I have a huge problem and I can not find a reason for it. I have two PCa, one in my house and the other in my workplace. The first one is an old notebook with Windows x86 and the other is Windows x64. I work in R in geospacial topics, and I have huge difference between the two system. The idea of my work is get the pixels of a raster (NetCDF) inside of a boundary of my whateshed (vector), so the code that I used is something like this:

 ncin <- nc_open(ncfname) #open the ncdf file
 #extract the variable 'precipitation' and latitud and longitud
 pp<-ncvar_get(ncin, varid = 'pr', start = c(1,1,1), count = c(-1,-1,1))
 lat <- ncin$dim$lat$vals   # latitude position
 lon <- ncin$dim$lon$vals  # longitude position
 # create a raster from ncin and transpose and flip the raster and extent
 r <- raster(pp, crs ="+init=epsg:4326" )
 rt <- t(r)
 rt <- flip(rt,direction=2)
 extent(rt)<-extent(c(range(lon), range(lat)))
 #open the shapefile or boundary (watershed), the shapefile is in CRS:4326
 myShape<-readOGR(dsn = myShapefileDir, layer = myShapefile)
 #crop the raster with the boundary
 crop = crop(x=rt,y = myShape, snap='in')

Then I have like the raster map is slipped in Windows x86 in comparison with Windows x64, you can look the difference in the next image (The left plot is created in Wx86 and the right in Wx64:

I know that right image is from GRASS GIS, but happens the same in R or GRASS GIS.

The correct image is the right, but I cant get the same results in Windows x86.

The netcdf file can download here at the end (it is call 'CR2MET_v1.4.2_pr_day_1979_2016') :

http://www.cr2.cl/datos-productos-grillados/?cp_cr2met=3

The boundary or shapefile can download from here (google drive, it is call 'polygon'):

https://drive.google.com/open?id=1N0w-mCYkPPk_kayTKnbkMJ8f3jsS4PuB


texts: Feature Layer will not refresh ArcGIS Javascript API

texts:
I am trying to integrate live point data (worker locations) from Workforce for ArcGIS into a Javascript web app that I am building.

The layer loads fine and displays on the map properly, it even refreshes as soon as I zoom in or out, but it does not refresh automatically no matter what I try. I have tried using .refreshInterval and when I log the FeatureLayer it has the property set correctly but no refresh happens.

  esriId.registerToken({
    server: 'https://www.arcgis.com/sharing/rest',
    token: sessionStorage.getItem('token'),
    userId: sessionStorage.getItem('user')
  });

  var portal = new Portal({
    authMode: 'immediate'
  });

  portal.load().then(function() {
    var query = {
      query: 'title:workers_ AND access:shared AND type:Feature Service'
    }
    portal.queryItems(query).then(function(queryResult) {
      workers = new FeatureLayer({
        title: 'Workers',
        refreshInterval: 0.1,
        portalItem: queryResult.results[0]
      });
      workers.when(function() {
        workers.makeTemplate();
        console.log(workers);
      });
      map.add(workers);
    });
  });

I have also tried using

setInterval(function() {
    workers.refresh();
    console.log('hi');
}, 6000);

This logs hi every time but still does not refresh the layer.
There are no errors of any kind in the console either. I have tried in both most recent Chrome and most recent Edge.

Does anyone know what the problem is/how to fix it?


texts: Is There A Term For The Old Maps That Had Lines Connecting Pin-Wheels?

texts:
Many old maps have lines like below. Is there a term for a map with that type of lines?

(from http://r-russell0912-dc2.blogspot.com/2011/01/pirate-maps-and-old-maps.html)

My best hypothesis is that the lines may be related to Rhumb lines and therefore Mercator projections.

Beautiful Example.


texts:
In arcpy.SearchCursor for row in rows I can return row.fieldname statically, but if fieldname is a variable, I don't know how to return the row.value of that variable dynamically.

for example:

    sortfield = something
    rows = arcpy.Searchcursor(feature)

    fieldlist = arcpy.ListFields(directory)
    for row in rows:

      for f in fieldlist:
        if f.name == something:
          id = row.something


texts: ArcGIS: Add XY Data isn't plotting in the correct location.

texts:
I've done this correctly so many times, I don't understand what's wrong.

The x and y columns I'm plotting are in decimal degrees.
Latitude is X and Longitude is Y.
I am plotting to Nad 1983 datum, which is consistent to my data frame and all my other layers. The points I'm plotting were taken in Nad 1983 datum.

They plot, but in the wrong location. What other things could be wrong?


texts: How to extract contours from .jp2 DEM dataset?

texts:
This will be a little long-winded but I'm including all information I consider relevant.

Problem:

I have received a data set from a municipality that is a 3-layer LiDAR JP2 DEM and need to extract ground contours at 1 m intervals for a specified area (more below).

Programs available:

QGIS 2.18 Las Palmas with LASTools add-on installed.
GRASSgis 7.4.4 with GDAL.
Civil3D 2017.
Land Desktop 2004.

Additional Information:

I normally work with LiDAR in .laz/.las formats and use LASTools through QGIS to extract the ground points with laszip and then create a .tin for use in Civil3D through Recap360 (.rcap import).

I can load the .JP2 DEM into GRASS with a proper spacial reference (UTM20), but it comes in as a 3-layer raster map with color naming being the only difference (R/G/B)

When I run r.contour to create a vector layer capable of being exported to .dxf, the returned vector polys are very noisy and unusable. I've tried a cheap work-around by increasing the minimum points required per contour but it's still unusable. I've also tried exporting the DEM .jp2 into a gtiff file but it's unrecognized by Civil3D (Civil3D takes .dem, .tif, .asc, .txt, and .adf for surface creation, but seemingly won't take a gtiff).

Assumed solutions:

I need to first merge the 3 band of scans, then extract the ground contours through vector .dxf, or
Combine the 3 bands and export as a singular .dem to be directly imported into Civil3D, or
Combine the 3 bands and export as a .shp, import .shp into LandDesktop, and then use the auto-created .tin file directly in Civil3D.

Any help is greatly appreciated. I can also upload the .zip of the scan to a cloud storage if someone wants to try their hand directly (it's public non-proprietary LiDAR set), but would like to know the solution regardless as I will likely need to do this again in the near future.


texts: Earth Engine: Downgrading image resolution

texts:
I am trying to downgrade the resolution of high-resolution precipitation product (DayMet) to low resolution using the solution posted here: https://gis.stackexchange.com/a/323143/33177

However, I am getting the following error:

Layer error: Image.reduceResolution: The input to reduceResolution does not have a valid default projection. Use reproject() first to set the base projection.

I have checked both the images have default projection.
Here is the code


texts: How can I access the array of predicted addresses from Google Maps once an address is entered?

texts: The Google Maps JavaScript Autocomplete API mentions that autocomplete returns an array of prediction objects.  How can I access this array without autocomplete?  I have a standard HTML form with address fields (like house number, street name, city, etc.).  The information is passed into a JavaScript function that concatenates the fields into a street address, then geocodes the address using the Google Geocoding API.


texts:
I've been reading about using Google Earth to query data (this article) and I'm interested in learning how to make one of these.

I just have a couple of questions:
1. Is this possible with the free version of Google Earth?
2. Can Google Earth connect to any database application (i.e. Access/SQL Serve)?
3. Is it a good query tool (in comparison to querying layers in ArcGIS or a relational database)?
4. How long does it take to make one of these tools?

Thanks for your help everyone!


texts: When I try to apply a downloaded GIS satellite (t=rimage to Quantum, it doesn't seem to open properly.  There must be a step-by step description on how this would be done effectively.  Can anyone tell me where to find it?


texts: How to know the date of adding specific layer to a map in ArcGIS?

texts:
Is there a way of knowing when specific layer have been added to MXD file?
I work with MXD file with three layer with the same name: "polygon":

i searched in the source tab (in the layer properties), but i didn't found this sort of data.


texts:
I have been given what I'm told are differential corrected gps points (DGPS) and asked to convert them to something usable in current GIS (e.g. Arcgis, Qgis).

They have names like R081320a.cor, R081323a.cor and assuming the file time stamps are accurate were created in 1999. Opening with a text editor reveals a binary format, with the first legible text as MCORR400 v4.55 W32, followed a little later with Automatically Created By GeoExplorer and Geo-Explorer Version 2.20.

Given the context, the GeoExplorer mentioned is probably the Trimble flavour, which needless to say I don't have.

What are my options for recovering these data?


texts: Map to KML not exporting all layers within data frame

texts: I have about 12 layers in my TOC originating from various sources (Folder, GDB and another folder). I would like to convert all this to a KMZ file. When I select "Layers" in the "Map to KML" tool data frame I expect that all of them should export. However only a few of them do. And whats even more strange: they originate from different sources... one is GDB and the other is folder. I would like to know how I can get the Map to KML tool to export everything all at once.


texts:
Please see my example OL3 map here of threatened ecosystems in South Africa: https://storage.googleapis.com/expanded-line-7365/DM%20selector%20%20with%20layer%20(demo)/index.html

Please use the search box to search for a District such as Eden, Ugu, West Coast, Sedibeng etc.

It is perfect how the search result highlights and zooms to the feature, but what I want to do is change the styling of this. For example I do not want the fill of the selected feature to become blurry, or possibly change the blue border to another colour? I have been searching through all the JS and CSS files but I have not been able to find where this information is set.

edit

The closest info I have found is this https://github.com/webgeodatavore/ol3-search-layer/blob/master/examples/control-ol3-search-layer.js

But when changing the stroke or fill values, nothing changes on my map.


texts: Querying a raster per pixel using a secondary raster as a condition

texts: So here's the problem.  I have two rasters containing bathymetry data: one containing the bathymetry of the sea floor, we'll call this raster1, and another containing the difference between two previous scans of the same area, we'll call this raster2.  I have to query raster2 such that any areas where the value is >10% of raster1's value will be highlighted.  I tried to do this through the Extract by Attribute in ArcGIS but could not figure out the SQL statement.  Any ideas on what the SQL query would look like?


texts: Reconnecting hyperlink in Dynamic Hyperlink mxd in ArcMap?

texts:
Is there a way to reconnect the Hyperlink in a Dynamic Hyperlink?

The situation is one of my workers was linking polygons to pdfs with Dynamic Hyperlink because each polygon had several pdfs that needed to be link.

The person sent me the mxd, and shapefile used for the Dynamic Hyperlink however the path to bring up the pdfs are not the same. The worker used his local drive and I trying to setup the link through our network so other workers can view the link pdfs.
Is there way to edit the path in a mxd file or the map properties to fix the link between the polygons and pdfs.

We are using basic ArcMap for the project.


texts: Getting the percentage area of intersecting polygons

texts:
I have the polygons for census tracts within New York City.
I also have the polygons for the parks.
What I would like is % of the census tract area, with an added 400 meter buffer, that is park.

I'm not sure if this is possible with Google Earth Engine. If it is, these are the steps that would be necessary.

Add a buffer to the census tract polygons
Calculate the area of each buffered polygon (A)
Calculate that area of park that intersects with each buffered polygon (Ap)
Divide Ap/A to calculate the percentage of park area within the buffered tract

Like I said, I'm not even sure if this is possible in GEE. Any pointers?

As requested, the data:

var parks = ee.FeatureCollection("users/tomlogan/NYC_parks"),
    tracts = ee.FeatureCollection("users/tomlogan/NYC_tracts");


texts: Adding layer with ArcObjects from Python script tool?

texts:
I would like to load an AFR file into a data frame in a currently open MXD and the IMap.AddLayer function is not working as expected if the data frame in which to load the AFR file is activate.
And I need to do all this with ArcObjects and from a Python script tool in a toolbox.

Here is what I do:

I apply a raster function template to a raster layer in a given data frame (a map).
I give the filename of the output raster dataset an ".afr" extension to ensure an in-memory temporary AFR file is generated.
I try to load the generated AFR file into my dataframe.

And it is Step 3 that is not working properly: the layer created from the .afr file quickly shows up and then disappears. Which means it seems to be added but does not remain in the data frame, in the table of contents.

At first I thought that the problem came from the .afr file I was not generated properly. But I further investigated and found out that even with a very simple layer addition the IMap.AddLayer function was not working if the targeted data frame (map) was active.

So I created a simple snippet that examplifies the problem so that people interested can reproduce the problem if they wish to.
Info: I am using Pierssen's Snippets module in order to access ArcObjects from Python. See Accessing ArcObjects from Python and Pierssen Snippet module for more information.

import Snippets
import comtypes.gen.esriArcMapUI as esriArcMapUI
import comtypes.gen.esriCarto as esriCarto

##############     MAIN     ##############

pApp = Snippets.GetCurrentApp()
pDoc = pApp.Document
pMxDoc = Snippets.CType(pDoc, esriArcMapUI.IMxDocument)

# Get list of data frames
pMaps = pMxDoc.Maps

# Get first data frame (= first map)
pMap = pMaps.Item(0)

filePath = r'C:\Temp\myRaster.tif'

pLayer = Snippets.NewObj(esriCarto.RasterLayer, esriCarto.IRasterLayer)
pLayer.CreateFromFilePath(filePath)
pLayer.Name = "myRasterLayer"
pMap.AddLayer(pLayer)

pMxDoc.ActiveView.Refresh()
pMxDoc.UpdateContents()

The raster layer is then supposed to be loaded in "DataFrame_1", which is the active data frame. But it is not.
 During the process it shows up very quickly but does not remain loaded in the data frame.
Also, I discovered that if "DataFrame_1" is not active, the layer is properly added. Weird!
Last important point: when this code is run directly from the Python window in ArcMap, it works properly whether the target data frame is active or not.

A few software details before ending: I'm using ArcGIS 10.3, with Python 2.7.

Did someone already encounter the same problem?


texts: Why ModelBuilder cannot see Summary Statistics field?

texts:
I'm having some problems adapting a few models inherited by me from a colleague.

Inside this model, we attempt to calculate Summary Statistics, but the model can't find the field that it needs to use. ModelBuilder refuses to show it to me, despite the fact that the Feature Class in this calculation possesses the field in question.

Has anyone seen weird problems like this one?

The original model was for 9.3 and I'm rebuilding it in 9.3.1. Any pitfalls I'm not aware?


texts: Polyline is not pointing to correct langitude.latitude data in Leaflet

texts:
I am trying to add polyline in my existing leaflet map. But unfortunately its not overlapping the correct coordinates.

this.testd = [
            [26.946345158021394, 45.30658810117468], [31.74737690870502, 58.27657611253403], [3.572373155835759, 50.87794848675696]
        ]

        this.pathLine = L.polyline(this.testd);
        this.allLayers.addLayer(this.pathLine);
        this.allLayers.addTo(this.map);

Please let me know how to make it work.


texts: Is there a Carto CSS Gallery which also contains code?

texts:
I am currently exploring the usage of CartoCSS and TileMill in one of our projects. I have seen the impressive gallery available at http://www.mapbox.com/tilemill/gallery/

The gallery looks impressive, but I am having trouble imaging how they achieved certain effects and styles. The Samples given in the Guide (for example here) seem too basic.

Are any detailed samples of the CartoCSS (rather than just the output) available somewhere?


texts: Identifying NoData pixels source in SoilGrids: current and legacy layers

texts:
I am trying to identify the origin of NoData pixels in the SoilGrids250 layers. i.e. I want to know which NoData pixels are inland water vs. sea water (the paper mentions using a sea mask) vs. permanent ice.

I am interested in both SoilGrids250 2017 (1.0) and SoilGrids250 2.0. This is because SoilGrids250 was used to make derived products for which I also need to identify inland water vs. sea water vs. permanent ice (e.g. Ross et al. 2018).

Are these layers available?

My understanding is that most of this processing (but not all) was performed using GlobeLand30 for SoilGrids 2017; but that website (http://www.globallandcover.com/User/Login.aspx) isn't functional for downloads (it's not possible to create a password + the forgot password page is down) and the emails of the data creators are invalid.

Note that I am specifically interested in identifying the pixels in the SoilGrids250 layers and not in improving upon these masks.


texts: Extending a segment in QGIS to check intersection with a polygon

texts:
I am trying to find the adress points (red dots) located in front of an obstacle (red polygon). The street is green and the pink line represents the street face (Si).

This is what I have so far, but another method would be acceptable as well.

Therefore I want to extend the blue line (which is the line from the adress point to the street face). And then check if the extension intersects with the polygon.

I am using the virtual layer function of QGIS but I can also simply do this with postgis. I tried to use this to write my query but I get mistakes (st_translate):

postgis, extrapolate a line

    SELECT
    adpt.id
    ,ST_TRANSLATE
        (adpt.id,
            sin(
                ST_AZIMUTH
                    (ST_ClosestPoint(si.geometry, adpt.geometry), adpt.geometry)
                )
            * (
                    ST_Distance(
                        ST_ClosestPoint(si.geometry, adpt.geometry)
                                )
                    * 10
                )
        ,
            cos(
                ST_AZIMUTH
                    (ST_ClosestPoint(si.geometry, adpt.geometry), adpt.geometry)
                )
            * (
                    ST_Distance(
                        ST_ClosestPoint(si.geometry, adpt.geometry)
                                )
                    * 10
                )
        )
FROM
    xl_adpt_ss_si ADPT
    JOIN xl_streetside si
        ON adpt.SI_ID = si.id
where
    ADPT.id='1136345'

;

Once I have my additional point, I want to use ST_MakeLine to draw the new line before computing the intersections.


texts: Add a field to output feature of a geoprocessing tool?

texts:
I am running the Watershed tool and RasterToPolygon tool in a loop for a number of points:

flowdir = "flow.tif"
flowacc = "flowacc.tif"
pnts = [arcpy.PointGeometry(arcpy.Point(x=0.0, y=0.0, ID = 12)), arcpy.PointGeometry(arcpy.Point(x=1.0, y=1.0, ID = 15))] #Create some test points
out = []

for i, p in enumerate(pnts):

    out.append("C:/temp/wspoly_%i.shp"% i)  #Create the output filename
    snap_pt = sa.SnapPourPoint(p, flowacc, tolerance = 100) # snap the point
    ws = sa.Watershed(flowdir, snap_pt) # calculate watershed raster
    arcpy.RasterToPolygon_conversion(ws, out[i]) # convert watershed raster to polygon

Ideally, I'd like the output to have an ID field taken from the ID attribute of each point.
Is there any way to do this?


texts: Patch processing of multiple HDF5 files using Python and Gdal (gdalwarp)

texts:
I am trying to extend this process

http://geoinformaticstutorial.blogspot.com/2014/06/reading-and-map-projecting-amsr-2-data.html

wherein a set of raster data with separate geolocation arrays are map re-projected using virtual rasters and gdal.

Goal: Create a global image mosaic from ~15 HDF5 "swath" grided data files. The final image must have rows = 4872 & columns = 11568. The latter is based on NSIDC EASE Grid Version 2 (http://nsidc.org/data/ease/ease_grid2.html)

The command line tags and options for gdalwarp are:

-geoloc -overwrite -co "COMPRESS=LZW" -tr 3000 3000 -srcnodata -9999 -dstnodata -9999 -r average-s_srs EPSG:4326 -t_srs '+proj=cea +x_0=0 +y_0=0 +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m +no_defs'

Each file has 6 data bands, thus at the end should get 6 different image mosaics. Files are ~2GB large.

Progress so far:

Able to open and read each file using h5py

indata = h5py.File(filename)
data_band_1 = indata["cell_sigma0_1"]
data_band_1 = np.where(data_band_1 <=0,-9999,data_band_1 )
data_band_1 [data_flg_band_1 == TRUE] = -9999

Save numpy arrays as Tiff using the following, but with out geotransform information

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
outFileName = mydates[iter_day]+'.tif'
print(outFileName)
gdalDriver = getGDALFormatFromExt(outFileName)
format="GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
newDataset = driver.Create(outFileName, inXSize, inYSize, numBands,gdal.GDT_Float32)
newDataset.SetProjection(srs.ExportToWkt())
newDataset.GetRasterBand(i+1).SetDescription(temp_name)
newDataset.GetRasterBand(i+1).WriteArray(temp_data)
newDataset.GetRasterBand(i+1).SetNoDataValue(nodat)
stat = newDataset.GetRasterBand(i+1).ComputeStatistics(0) newDataset.GetRasterBand(i+1).SetStatistics(stat[0], stat[1], stat[2], stat[3])
newDataset.FlushCache()
newDataset = None

Created virtual rasters (.vrt) similar to the method in the link above. But replaced the source data with the new Tiff.

Question: How can I merge all the .vrt files to a raster with a specific size (e.g., known rows and columns)

Note: the lat/lon steps are not necessary uniform (can't setup a world file)


texts: Extract Differential Interferometry from SENTINEL-1 image

texts:
I downloaded some SLC images from scihub.esa.int . And I want to extract Differencial Interferometry from them by using SNAP-ESA software.

How can I do it ?


texts: How to get max value from 3 vector files?

texts: I'd like to know how to the get max value from 3 vector files. Originally, the vector file are interpolated images. (First, point data are interpolated into surface data by IDW method and then exported into vector file). Now I need a vector file with max value obtained from the 3 vector files. When I work on raster image to get max value, firstly, I always do layer stack operation on all images from which max value will extract. After that, I get max value image from stacked image. So what I want to ask is do I need to do layer stack on 3 vector images? If so, what is the operation name/tool in ArcGIS and how can I find max value from stacked vector files?


texts:
So I have a layer with thousands of gaps and overlaps that I want to tidy up, but when I go to run v.clean in order to do so, it just sits there at 0%. I have no warnings or errors in the system, and QGIS itself is not responding. Normally when this happens, I just leave it alone and it sorts itself out. However it's been over two hours now and nothing has changed.
I'm running the latest version of QGIS and GRASS, and there's plenty of unused RAM and disk space the program could be using, but it's not. What is going on and how can I fix it?


texts: Mapping multiple categories that overlap, but cannot get layer to appear on top consistently ArcMap 10.5

texts: I am trying to map overlapping practices using unique values in ArcMap 10.5. I have a practice for planting and practice for site prep. Planting is a solid fill blue but for some reason the prep practice that is set to crosshatch 10% is only appearing on 1/5 of the polygons that it should. I am a novice so I will attach a picture if someone can help. The area with blue should all be showing cross-hatch. I tried adjusting the advanced symbol levels, but that just got rid of all of the crosshatch. I am not sure exactly how the symbol levels work.